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Abstract 

To model the radiative evolution of extreme mass-ratio binary inspirals (a key target of the LISA mission), 
the community needs efficient methods for computation of the gravitational self-force (SF) on the Kerr 
spacetime. Here we further develop a practical 'm-mode regularization' scheme for SF calculations, and 
give details of a first implementation. The key steps in the method are (i) removal of a singular part of the 
perturbation field with a suitable 'puncture' to leave a sufficiently regular residual within a finite worldtube 
surrounding the particle's worldline, (ii) decomposition in azimuthal (m-)modes, (iii) numerical evolution of 
the m-niodcs in 2+ ID with a finite difference scheme, and (iv) reconstruction of the SF from the mode sum. 
The method relies on a judicious choice of puncture, based on the Detweiler-Whiting decomposition. We 
give a working definition for the 'order' of the puncture, and show how it determines the convergence rate 
of the m-modc sum. The dissipative piece of the SF displays an exponentially convergent mode sum, while 
the TO-modc sum for the conservative piece converges with a power law. In the latter case the individual 
modal contributions fall off at large m as m~" for even n and as m~"+^ for odd n, where n is the puncture 
order. We describe an m-mode implementation with a 4th-order puncture to compute the scalar-field SF 
along circular geodesies on Schwarzschild. In a forthcoming companion paper we extend the calculation to 
the Kerr spacetime. 
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I. INTRODUCTION 



The closely-related notions of "self force" [U Ej , "radiation reaction" j3] , and "radiation damp- 
ing" [3] have a long and interesting history in physics. A motivating example arises in classical 
electromagnetism when one considers a point particle undergoing acceleration. An accelerated 
charge produces electromagnetic radiation; hence the particle loses energy and must therefore ex- 
perience a braking force. The braking force may be interpreted as arising from the interaction of 
the particle with its own radiative field. This interpretation is not straightforward mathematically, 
since the electromagnetic field is formally infinite at the particle. Dirac [3] showed how to remove 
a divergent (time-symmetric) component of the field, to isolate the finite (non-symmetric) part 
responsible for radiation reaction. 

The self- force idea now finds a modern application in the study of Extreme Mass Ratio Inspirals 
(EMRIs), which are a key target of the LISA mission [SH?]. An EMRI is a special case of the 
gravitational two-body problem, in which a compact body (e.g., a stellar-mass black hole or neutron 
star) of mass fi is gravitationally bound to a massive black hole of mass M, such that the mass 
ratio fi/M is very small. If n/M is vanishingly small, the smaller body follows a geodesic in 
the spacetime of the larger body [8|. For a small but non-zero /x, the smaller body perturbs 
the spacetime geometry at 0{p/M), and the resulting back-reaction force is referred to as the 
first-order "gravitational self- force". 

It was appreciated long ago [H |9l [10] in the electromagnetic context that self- force calculations 
on curved spacetimes are more challenging than in flat space. In principle, on a curved spacetime, 
the self force (henceforth SF) depends on the entire past history of the motion. Hence the challenge 
is not merely to derive formal expressions for the SF, but also to implement practical schemes for 
their evaluation (see [1112] for reviews). 

Following the foundational work of Dirac [3J and DeWitt and Brehme [9] on the electromag- 
netic SF, a key step forward came in 1997 with the derivation of equations for the (first-order) 
gravitational SF |1H I12j. now known as the "MiSaTaQuWa" equations. Alternative derivations 
and extensions have appeared over the subsequent years, for example in the works of Detweiler and 
Whiting [13], Gralla and Wald [8], Harte [HlIS], Gralla et al. [IS], and Pound [I?]. Equations for 
the scalar-field SF [U [18] and the electromagnetic SF [9] [lOl llSl [1^ have also been obtained. 

It is a non-trivial task to compute the gravitational SF from the MiSaTaQuWa equations, and 
first they must be cast into a form amenable to practical computation. One standard method is 
the so-called "/-mode regularization scheme" , outlined in |20[ [2T] . This scheme has been applied in 
a range of studies of SF on the Schwarzschild spacetime, for example, the scalar SF for radial infall 
[22] , circular orbits [23l - [26] and eccentric orbits [271 HH] ; electromagnetic SF for eccentric orbits 
|29j : and the gravitational SF for radial infall |30j . circular |3H - t35j and eccentric orbits |36j . Other 
approaches under development include [37H32] . 

There are promising signs that the SF programme is approaching maturity. For instance, 
gravitational SF results are now being compared against results of post-Newtonian theory |43H47j , 
and used to calibrate functions in the Effective One-Body (FOB) theory in the strong-field regime 
\45\ I46j . The shift in the innermost stable circular orbit (ISCO) due to the conservative part of 
the gravitational SF was recently computed [58]. SF results are also being used to inform data 
analysis strategies for the mock LISA data challenge [l9]; an improved understanding of strong- 
field SF-related phenomena such as resonances |50l |5T] will be undoubtedly prove central to this 
effort. There now arises the possibility that SF results will soon be meaningfully compared against 
Numerical Relativity simulations, which are pushing into the intermediate mass-ratio regime [521 - 
[55]. 

Whilst most studies thus far have assumed the central black hole to be non-rotating 
(Schwarzschild- type), it is reasonable to expect that in astrophysically-relevant EMRIs it will 
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be rotating (Kerr-type). At the time of writing, few calculations have been attempted for the 
more technically-demanding Kerr case. An exception is the recent work in Ref. [^H]) in which the 
/-mode regularization scheme is applied in the frequency domain to compute the scalar SF for 
equatorial circular orbits on Kerr. Unfortunately, the lack of separability of the gravitational field 
equations means the l-mode scheme cannot be applied in a straightforward way to gravitational 
SF calculations on Kerr. This motivates the development of alternative methods. 

SF calculations may be performed in either the frequency or time domains. In the frequency 
domain, the SF is reconstructed from a sum over frequency modes, with a spectrum of frequencies 
which are integer multiples of the fundamental orbital frequencies. The frequency-domain approach 
works well for circular orbits and low-to-moderately eccentric orbits (e < 0.7 [57J). It can give 
highly accurate results, because a complete decomposition (into frequency and angular modes) 
leaves one with an ordinary differential equation for the radial functions, which may be solved to 
high precision. 

The frequency-domain approach has limitations, however. It is not suitable if the field equations 
cannot be separated (as in the important case of metric perturbations on Kerr in Lorenz gauge), 
or if the orbit has high eccentricity [58], is highly generic (on Kerr), or is unbound. Importantly, 
frequency domain methods seem much less well suited to the challenge of evolving an orbit per- 
turbed by a SF in a self-consistent manner. This has motivated the development of a range of 
time-domain approaches. 

The 'm-mode regularization method' introduced in [591 [60] provides a general framework for 
time-domain SF calculations in axisymmetric spacetimes (such as Kerr), which may be applied 
in the scalar, electromagnetic and gravitational cases. Let us briefly examine the similarities and 
differences between the I- and m-mode approaches, which are both based on a decomposition in 
angular modes. Whereas the l-mode decomposition in spherical harmonics results in 1 + ID modal 
equations, the m-mode decomposition in azimuthal modes leads to 2 + ID equations. Here lies 
a key difference: field modes in 2 -|- ID are divergent on the worldline, whereas 1 + ID modes 
are continuous (albeit not differentiable) there. This obviously poses a challenge to numerical 
schemes, and motivates 'regularization' of the m modes with an analytically-determined 'puncture 
field' which is used to remove the divergence. In Ref. [59] it was shown that, with a simple ('1st- 
order'; see below) puncture field, the 2 + ID field modes could be evolved numerically. The original 
idea behind |i59j was to use the m-modes to obtain (via integration over 9) the /-modes which 
are needed as input to the standard /-mode regularization scheme. However, it was later shown 
in [60] that, with an improved ('2nd-order') puncture field, the SF could be recovered directly 
from a sum of the gradients of the m modes themselves. The current work describes the first 
implementation of this idea: we use the m-mode regularization scheme to compute the scalar-field 
SF in Schwarzschild. In a companion paper we will describe a scalar-field implementation on the 
Kerr spacetime. 

The aim of this work is to lay the necessary foundations for, and to demonstrate the com- 
putational feasibility of, accurate m-mode 2 + ID time-domain SF calculations using high-order 
punctures. The punctures employed in [59l|60] are motivated by the Detweiler- Whiting split [13j of 
the retarded field into 'radiative' (R) and 'singular' (S) parts. We give a method for constructing 
puncture fields from finite-order local expansions of the S field (see also Ref. [61] ) , and demonstrate 
that the order of the expansion directly affects the convergence rate of the m-mode reconstruction 
of the SF. We carefully describe the features of the m-mode scheme that will be common to all 
future implementations, such as the puncture formulation, the puncture order, the m-mode conver- 
gence rate for dissipative and conservative parts, the worldtube formulation, initial and boundary 
conditions, and the modelling and mitigation of various sources of numerical error. In this work 
the Schwarzschild spacetime serves as a testing ground in which to explore the features that do 
not depend on the precise form of the field equations. In the companion paper, we describe an 
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implementation that explores the issues specific to Kerr, such as stability of finite-differencing 
schemes. 

Other time-domain approaches are under active development. Vega et al. [41j have outlined a 
framework for time-domain SF calculations in 3-l-lD. They have demonstrated that, with a little 
modification, codes originally written for applications in Numerical Relativity can be used for SF 
calculations. Furthermore, they have computed the scalar SF for circular orbits on Schwarzschild to 
within ~ 1% accuracy. Pushing the accuracy towards the one-part-per- million benchmark achieved 
by Thornburg |62] in a I+ID time-domain code with adaptive mesh refinement represents a consid- 
erable challenge, and an ongoing community effort is underway. We believe the m-mode approach 
represents a competitive alternative to the 3-l-lD scheme, since it exploits the axisymmetry to 
achieve substantial gains in computational efficiency. Of course, the price to pay for decomposition 
is that the SF must then be reconstructed from a sum over modes; but, as we aim to show here, 
the convergence of the mode sum is now well understood. 

The paper is arranged as follows. In Sec.|TI]we outline the theoretical basis of our approach. Here 
we cover the Detweiler- Whiting split into 'radiative' and 'singular' fields (II A) which motivates 
the 'puncture' scheme (II C), the definition for 'puncture order' (IID) and its effect on mode sum 
convergence properties (II G), and the worldtube formulation (II H). In Sec. Ill we give details 
of the first 2-l-lD implementation for a SF calculation, for circular orbits in Schwarzschild. We 
describe the calculation of the puncture function and the effective source at orders 2, 3, and 4 



(IIIB), the code architecture, the finite difference scheme and numerical stability (|IIIC), and the 

In 



method for reconstructing the SF from data extracted from multiple 2+lD 'runs' (HID 



HIE). 



Sec. IV we present a sample of numerical results. After exploring the various sources of error and 
strategies for error mitigation (IV A TVB[ ), and numerically testing various predictions of Sec. |ll| 
we present in Sec. IV D the m-mode SF results, which we compare against the frequency domain 
results of Ref. 



We conclude in Sec. |V] by reviewing progress and outlining themes for future 
work. The Appendices make explicit some of the more elaborate technical parts of our calculations. 
Throughout we adopt the metric signature (—1, 1, 1, 1) and set G = c = 1. 



II. THEORETICAL EXPOSITION 

Consider a test particle with scalar charge q moving along a worldline 7 in the vacuum exterior of 
a black hole. Neglecting SF effects, we assume that the worldline 7 is a geodesic on the background, 
parameterised by z{t) = z^{t) where r is the proper time, with a tangent vector = dz^/dr. The 
charge acts as a source for a scalar field which satisfies the minimally-coupled Klein-Gordon 

equation, 

= V^V^$ = S{x), (1) 

where denotes the covariant derivative with respect to the background metric g^u- The source 
term is given by 

/oo 
5^ {x - z{t')) dr', (2) 
-00 

where p is the charge density, g is the metric determinant, and S'^ is the four-dimensional Dirac 
delta distribution. The retarded solution to this wave equation may be expressed as 

/poo 
G,et(x, x')pix')d^x' = q / G,et(x, z(r')) i-giz))-^/^ dr', (3) 
J —00 
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where G^eti^jx') is the retarded Green function defined by 

a^Gretix, X') = -4^5^ {x - x') , (4) 

with appropriate retarded boundary conditions. The retarded Green function has the Hadamard 
form |63] 

Gret(x,x') = e_(x,x') [U {x , x')6{a) + V {x , x')Qi-a)] , (5) 

where U and V are regular symmetric biscalars, and a = a{x,x') is Synge's world function [64j . 
defined to be half the squared geodesic distance between spacetime points x and x' , with a being 
negative (positive) when x and x' are connected by a timelike (spacelike) geodesic, and zero in the 
null case. Here 0(— cr) is the Heaviside step function, and 0_(x,x') is unity if x' is in the causal 
past of X and zero otherwise. Note that here we have adopted the sign convention for V{x) of [l], 
which is opposite to, e.g., 165], 

Note that the Hadamard form, Eq. ([5]), is only valid when x' may be connected to x by a unique 
non-spacelike geodesic. More precisely, Eq. ^ is valid if and only if x and x' lie within a convex 
normal neighbourhood [66]. Particularly, in the presence of a black hole, this condition is not valid 
for points x' on the worldline in the 'distant past' of x (see, e.g., |38] for a discussion) and so Eq. ^ 
is of restricted utility. 

The particle obeys the equation of motion 

n,V^(/xnn = Fi^eif, (6) 

where -Fggjf is the scalar SF. The component of F^^^^ orthogonal to gives rise to a self-acceleration; 
the component tangential gives rise to a change of mass (in the case of circular motion it turns 
out there is no tangential component and thus no change of mass). From a naive application 
of a Lagrangian principle [67|, one would expect the scalar SF to be obtained as the gradient of 
the retarded scalar field, gV^$ret- However, this expression of course becomes meaningless when 
evaluated on the worldline, where $ret and V^<I>ret diverge. A more careful and considered analysis 
is needed. It has been shown [U |18] that the scalar SF along a geodesic on the vacuum exterior of 
a Kerr black hole is given by the integral 

F-if(T) = ^li^m J^^^ V^G^et {x, ^(r')) |_,(,,) dr' . (7) 

It is difficult to evaluate this expression, because it involves an integral over the entire past history 
of the particle's motion. Next we consider an alternative expression more amenable to practical 
computation. 



A. Detweiler- Whiting R S decomposition 

In an influential work on classical electromagnetism in flat spacetime, Dirac [3] showed that 
the physical radiation reaction force could be obtained if one removes a certain singular and time- 
symmetric component from the physical (retarded) vector potential, to leave behind a 'radiative' 
field. Detweiler and Whiting |13j elegantly extended Dirac's argument to curved spacetimes (for 
the scalar, electromagnetic and gravitational cases). In particular, they have shown that the scalar 
SF equation ([T]) is alternatively obtained by taking the derivative on the worldline of a certain 
radiative/regular (R) field, 

Ff{T) = q lim V^a>^(x), (8) 

X—>-Z(T) 
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where ^r{x) is an homogeneous solution of Eq. ([T]) (i.e. with 5 = 0) given by 



(9) 



Here the symmetric/singular (S) field is a particular solution to the sourced wave equation given 
by 



Gsix,z{T))dT, 



where the symmetric Green function Gs is defined through its Hadamard form 

Gsix,x') = ^ [Uix,x')6ia) - V{x,x')Qia)] . 



(10) 



(11) 



Here U and V are the same biscalars that feature in Eq. (5|. By construction, Gs is zero inside 
the future and past light-cones. Inserting Eq. (11) into (jlO) leads to 



Uix,z{T)) 

2a 



+ 



Tret (a;) 



U{x,z{t)) 



2d 



TadvCa;') 



"^adv 



V{x,z{T))dT, 



(12) 



where 2;(Tret) and z{ts,^^) are the points on the worldline in the causal past and future of x that 
are connected to x by a null geodesic, and Tret and Tadv are the corresponding proper times along 
the worldline. Note that Ti.et(a^) and Tadv (a;) are non-smooth functions of the field point x that are 



not differentiable on the worldline, and that strictly speaking (12) is only well-defined if x and z 
lie inside a convex normal neighbourhood. 



B. Dissipative and conservative parts of SF 

The SF can be split into 'dissipative' and 'conservative' parts, as follows. First, let us introduce 
the 'advanced' field 'I'adv; which is defined in an analogous way to $ret; i-e., via Eq. ([s]) with a Green 
function Gadv(2;, x') governed by Eq. ^ with appropriate advanced boundary conditions. Next we 
may define retarded and advanced radiative (R) fields via $5?* = *^ret — and <I>^^ = <l*adv — ^s- 
Note that the same singular (S) field is used in both cases, since it represents the 'symmetric' part 
of the field. Then we may define the 'conservative' and 'dissipative' combinations 

^cons ^ 1 (^^ret ^ ^adv^ ^ 1 (^^^^ ^ ^^^^ _ 3$^^ ^ (^3) 
^diss ^ 1 ^^ret _ ^adv^ ^ 1 (^^^^ _ ^ (^4) 

and the corresponding 'conservative' and 'dissipative' parts of the scalar SF follow from 

^cons/diss ^ ^ lim v^$""'^/<^'^^ (15) 

x— 5-z(r) 

A key point here is that the dissipative component is found from the field difference ^rot — ^adv, 
which is known to be a smooth function even on the worldline. In other words, provided we 
can obtain the advanced field the dissipative component of the SF does not require regularization. 
On the other hand, to compute the conservative component one requires knowledge of the S field, 
which is singular on the worldline. 
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C. Puncture scheme 



Computing <I>5 is not straightforward. Unfortunately, a closed form expression for the biscalars 
U and V is not available, and generally one falls back on approximation methods. Fortunately, 
methods exist to expand the biscalars U and V as covariant series in a, and ultimately as series 
in coordinate separations. The Hadamard-expansion method is now well advanced for several 
spacetimes of physical relevance, such as Schwarzschild and Kerr [65 1 l68l - [70] . 

Knowledge of the Hadamard expansion of $5 allows the introduction of a 'puncture' {V) field 
^■p. To be of use, must be defined 'globally' (or at least everywhere within a region surrounding 
the worldline), and in the vicinity of the worldline it must have the same local behaviour as the 



S field of Detweiler and Whiting [Eq. (12)], up to a certain order (to be made precise below). We 
then define a 'residual' (7^) field to be 

^n{x) = <^>rct{x)-<^v{x). (16) 

The residual field $7^ will have the same local expansion as ^ji, up to a certain order, and (if the 
order is sufficient) the scalar SF may be obtained from evaluating the derivatives of ^-ji on the 



worldline (see Sec.|IIF). This basic idea is also applicable in the electromagnetic and gravitational 



cases [60j; in this work we focus on the scalar-field case. 



D. Order of the puncture function 

Let us now give a working definition for the 'order' of the puncture function, with reference 
to local coordinate expansions. Our classification of order should in fact be independent of the 
choice of coordinates, provided one works with a sufficiently regular coordinate system. A similar 
definition of order from a covariant point of view is given in [61J. 

We begin by making a rather subtle distinction between a field ^-p{x; z) which is defined locally 
in the vicinity of a particular point z on the worldline, and a puncture function <^p(a;; 7) which must 
be defined everywhere in the vicinity of the worldline 7. Likewise, we should distinguish between 



the global residual field, defined by Eq. ( 16 ), and a local residual field ^-j^ix; z) = (^j-ct{x) — ^p{x; z) 
which is again defined with reference to a particular point z on the worldline. In Sec. [HE] we give a 
practical scheme for promoting a 'local' expansion <I>^ to a 'global' puncture field <^-p in the simple 
case of circular orbits. 

Denote a field point by x'^, denote a worldline point by z^, and define the coordinate differences 
5x^ = x^^ — z^. It is convenient when discussing the order of the puncture to introduce a scaling 
parameter A and new coordinates 5x^, through 

5x^' = X5x^'. (17) 

Taking the limit A — >■ with fixed 5x^^ is equivalent to approaching the point from a specific 
direction. 

A lowest-order puncture was given in |59| : 



V 

and 



$gJ(A5x) = (?/€[!], where e^^^ = \X\Sl'\ (18) 



So= {g,,, + u^u,)\Jx^'5x\ (19) 
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Subtracting ( 18 ) from the retarded field leaves a residual field <I>^^ which is C ^ in the sense that 
it is bounded but discontinuous as 6x — )■ 0, viz. 

lim $W(A(5x)/ lim <^^l\x6x) (20) 

in general. 

The next-order puncture was obtained in [60J, and may be written 

^^'='?/e[2], where epj = |A| (5o + A^i)^/^^ (21) 

and 

Si = {uxu.T^^, + 5/..,a/2) \jx^'5x''5x'', (22) 

where T^^^u are Christoffel symbols for the background metric. An alternative puncture of the 
same order may be defined as 

so that <|)^''^^*^ — = 0{\X\). Subtracting either puncture from <I>ret leaves a residual, i.e., a 

[21 

function ^ ~ which is continuous but not differentiable at 6x = 0. 
n 



Starting with Eq. (12), the Detweiler- Whiting S field may be expanded as 

where Fk{5x) are polynomials in 5x of order 3/c. We will call <1>^' an "nth order puncture" if 

-^s = (|A|A"-2) and lim U^;^ - $5) / (|A|A"-2) ^ ^25) 



It follows that 

=^fi + 0(|A|A"-2). (26) 
Since is a smooth function, the residual field is C"~^. Hence, for example, a 2nd-order 

[21 [31 

residual <I>i- is continuous but not differentiable, and a 3rd-order residual ^ ~ is both continuous 
and differentiable. 



E. Global definition for the puncture function 

The covariant expansion method developed by Ottewill and Wardell [521 HO] enables one to 
compute expressions for nth order 'local' punctures expressed in terms of coordinate differences, 
i.e. <I>p^(5x'^). In this paper, we discuss 2nd, 3rd and 4th-order punctures; in principle, higher 
orders are possible (see Ref. [61] for a discussion). 

Given a field point x, we are free to choose the worldline point z to lie anywhere on the worldline 
between zir^etix)) and z{Tg^A^{x)). In order to obtain a puncture function which is globally-defined 
(or at least defined within the vicinity of the worldline for all t), we must allow z to become a 
function of x. 
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There is more than one way to relate x to z. For instance, we could choose z to be the point 
on the worldline that is connected to x by a spacelike geodesic orthogonal to the worldline at the 
point of intersection. Although a natural definition, this makes z a rather complicated function of 
X. A simpler approach (and that used in e.g. [59]) is to set the coordinate time of z to be equal to 
the coordinate time of x, z^ = = t. This is a coordinate-dependent construction; henceforth we 
will work in the Boyer-Lindquist system {t,r,9,(p}. Then 

dt = 0, 6r = r-rp{t), 6e = 9-ep{t), 5cp = if - ipp{t), (27) 

where rp{t),9p{t),ifp{t) are coordinate functions describing the worldline, and a globally-defined 
puncture function for use in our scheme is 

$^](xn = $^^(0,r - rp{t),e- Opit),^- ^p{t)). (28) 

A puncture function of the form (|28|) will generally be ill-behaved at spatial infinity, and possibly 



elsewhere; we deal with this problem in Sec. II H There is still arbitrariness in the definition of 
the puncture function since the only requirement is that the puncture field has the correct 



expansion [see Eq. (25)] in the vicinity of the worldline. For example, within the m-mode scheme 
we are motivated to replace 5ip with a smooth periodic function f{Sip). This would not change the 
order of the puncture if /((/?) = Sip + 0{(p^^^). 

F. Self force from the residual field 



Let us define the residual field through Eq. (16) with a 'global' puncture field (28), and now 



consider its gradient near the worldline, i.e. the quantity 

V^cdN =v^<I>,j + 0(|A|A"-3). (29) 

It is clear that this quantity is only guaranteed to be well-defined on the worldline if n > 3, i.e. if 
we use a 3rd-order puncture or higher. It is also clear that in this case, the gradient evaluated on 
the worldline leads to the correct SF via Eq. ([s]), i.e. 

Ff'{r)=q hm V^cI>S^^^l(x). (30) 

However, this is not the complete story. If we employ a 2nd-order puncture, the gradient is 
discontinuous at the worldline, 0(|A|/A). In other words, it depends on the direction in which the 
worldline is approached. Nevertheless, it was shown in |60| that the SF constructed from a sum 
over azimuthal m modes is in fact well-defined and correct. A related fact is that the convergence 
rate of the m-mode sum depends in a particular way on the order of the puncture function, as we 
now begin to discuss. 

G. TO-mode decomposition and mode sum convergence 

We may take advantage of the azimuthal symmetry of the background (i.e. the Kerr spacetime) 
to decompose into m-modes, i.e. 

oo 

Q{t,r,9,v) = Yl Q"'it,r,e)e'^^, (31) 

m=—oo 

Q"'{t,r,9) = ^ Q{t,r,0,ip)e-'"''^dip, (32) 
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where Q may be any member of the set {^>ret; ^p^, 'S'eff} [here is the effective source to 
be defined in Eq. ( |43| ) below]. For later convenience, we define the 'total m-mode contribution' 
Q'^{t,r,6,(p) (or 'modal contribution' for brevity) to be the real quantity given by 

^ ^*''^'^''^^- \ Q™(t,r,0), m = • 

The value of the radiative field ^r{z) at a point z on the worldline may be found from a sum over 
the modal contributions of the residual field, 

^n{z) = Y: ^^n-'^-'iz). (34) 

m>0 

The SF is reconstructed from the gradient of the residual field modes evaluated at the worldline 
point z, 

Ff\z) = Y,F7i^). where F™(z) = g lim V^l>fe>^l™. (35) 

m>0 



1. Exponential convergence of the dissipative SF 



In Sec. II B we described the split of the SF into dissipative and conservative pieces. The 
dissipative piece F^^^^^ is found from the gradient of ^'^^^^ , where ^'^^'^^ is formed from the difference 



between retarded and advanced fields [see Eq. (14)]. The difference ^'^^^^ is a smooth (C°°) function, 
hence its m-mode contributions decay faster than (where k is any positive integer) in the limit 
m — )• oo. We will call this behaviour 'exponential convergence'. However, in order to construct the 
gradient of ^'^^^^ in practice we need a method for calculating the gradient of the advanced field 
explicitly. In the case of eccentric geodesic orbits in the equatorial plane, it is straightforward to 
obtain V^<I>adv on the worldline by making use of the convenient symmetry relation [see Eq. (2.80) 
in [50]] 

V^^>adv(r,Ur) = -efj,Vf,^ret{r,-Ur), (36) 

where = (1, —1, —1, 1) and there is no summation over /x. In other words, by identifying a point 
on the orbit conjugate to the point of interest (i.e. a point at the same radius with equal and 
opposite radial velocity Ur), we may obtain the gradient of the advanced field directly from the 
gradient of the retarded field (see, e.g.. Sec. HE in [36j). 

Let us consider the case of circular orbits in the equatorial plane in a little more detail. In this 



case, since Ur = 0, it follows via (36) that V^<I>adv = iV^<I>ret) where the gradients are evaluated at 



the same point on the worldline. Here the plus sign corresponds to the r component, and the minus 



sign to the t and (/? components. It follows immediately via (13) and (14) that the conservative 
part of the SF is purely radial, whereas the dissipative part has components only in the t and <f 
directions, i.e. F^^^^ = (^FI^^'^^,F^°^^,0,F^^^^). Now, let us consider the (dissipative) t component of 
the SF; by the above discussion it follows that 

F[^ = q lim Vj4>5>^''" = q lim V*!'^^^^''", (37) 

where $diss,m ^^iq m-mode contribution to Since the field ^'^^^^ is a smooth function, 

it follows that the quantity on the right-hand side converges exponentially fast with m. Hence 
the SF modes F[^ must also convergence exponentially fast with m. A similar argument follows 
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immediately for the (dissipative) (/9 component of SF, which in fact is related to the t component 
via 



pself 



(38) 



where uj is the frequency of the circular orbit. Note that the argument for exponential convergence 
is independent of puncture order, and so the modal contributions -F™ and F"^ are also independent 



of the puncture order (we check this in the numerical implementation of Sec. IV, cf. Fig. 15). 



2. Power-law convergence of the conservative SF 

Let us now investigate the m-mode convergence properties of the conservative part of the SF. 
A careful analysis of convergence for the 2nd-order puncture scheme was given in [60]. In this 
section, we eschew a formal analysis in favour of a heuristic analysis which illustrates the key 
features. It remains for us to demonstrate that these features are supported by the results of our 



specific implementation, which we do in Sec. IV 



Let H^^'{(p) be a smooth function on — vr < 93 < vr, and across ip = — vr, vr (i.e. all its derivatives 
match there), except at 99 = 0, where it admits the local expansion 

^"■'-^|(^/. (39) 

Here are constant coefficients and hn,hn+i 7^ 0. The function H^^^ip) is akin to the nth-order 
residual field ^>^' . Note '"l is continuously differentiable n — 2 times everywhere, but its (n — l)th 
derivative has a jump discontinuity at 99 = 0. ff'"! has a mode-sum reconstruction of the form 



00 

m>0 



The m-mode contribution ^N'" [defined as in Eq. (33)] is shown in Appendix |a] to have the 
following asymptotic behaviour in the limit of large m: 

fVMm 2/i„ j i-)"/^ COS mip, neven, 
vrm" I (— ) 2 sinmip, n odd, 

+ ^ X I (-rj^-nm^, n even, ^ ^ / 
7rm"+^ 1 ( — ) 2 cosrrup, n odd, V / 



In the large-m limit, the m-mode contributions (41) decay as ~ 1/m", in general. However, 
at the irregular point ip = (i.e. on the worldline) they decay as ~ l/m""*"^ if n is odd. In other 
words, for odd n the mode-sum reconstruction is "one order more convergent" than would be 
naively expected. 

The toy model illustrates a well-known feature of Fourier theory: the smoother a function, 
the more rapid the convergence of its Fourier series. It also demonstrates a less obvious feature: 
odd-order punctures (i.e. n = 1, 3 etc.) will generate a residual field whose m-mode contributions 
decay one order faster than expected, i.e. ~ l/m""^^. 

Now let us consider the gradient of the residual field (giving the SF), which is in general one 
order less differentiable than the field itself. Odd-order punctures give modal contributions to 
V^<I>^' which decay as expected, i.e. as ~ l/m"~^. Even-order punctures (i.e. 2nd, 4th, etc.) give 
modal contributions for V^<I>^^ which decay one order faster than naively expected, i.e. as ~ l/m"'. 
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Puncture Order 




smoothness of 


S [rijm 


V 


^cff 


n — 1 
n = 2 
n = 3 
n = 4 


|A|/A 
|A| 
|A|A 
|A|A2 




_? 

TO 

—4 

TO 

—4 

TO 


_2 
TO ^ 

— 2 

TO 

—4 

TO 


|A|/A^ 
|A|/A2 
|A|/A 
|A| 



TABLE I: Differentiability and large-TO behaviour of key quantities. As discussed in the text, the convergence 

[nl 

rate of the TO-mode series is related to the smoothness of on the worldline, which depends on the order 
of puncture used. The parity of the puncture order is important; convergence of the modal contributions 
improves in jumps of to~^ every second time the order (n) is increased by one. The 4th and 5th columns 
show the large-m power-law behaviour of the individual to modes, for the residual field and the conservative 
part of the SF, respectively. The modes of the dissipative part of SF converge exponentially fast (not shown; 



see text). The final column shows the smoothness of the effective source S^s (to be defined in Sec. II H) 
near the worldline. 

The expected convergence behaviour of key quantities in our problem is summarised in Table |I} 
Let us note in passing that to obtain a mode sum for the conservative part of the SF whose terms 
fall off as m~^, we must use a 4th-order puncture; the 3rd-order puncture gives only modal 



fall off, just like the 2nd-order puncture. In Sec. IV B we give numerical evidence in support of 



this assertion (cf. Fig. 14). 



H. Effective source, worldtube formulation and modal equations 

The residual field is governed by an inhomogeneous wave equation, 

□cI>M=5S, (42) 

fnl 

where the effective source S^^ is given by 

s'§{x)^S-U^f. (43) 

The behaviour of near 7 depends on the order (n) of the puncture field. The puncture field is 
such that the distributional component of the original source (i.e. the delta function) is eliminated. 
In the vicinity of the worldline, (expressed in terms of coordinate differences = X6x^) has a 
local expansion starting at OdAlA""^). In other words, is divergent for n = 1, 2, discontinuous 

for n = 3 and continuous (C^) for n = 4. An illustration of a typical 4th-order effective source 
close to the worldline is shown in Fig. [Tl for a particle on a circular geodesic orbit. 

[nl 

The effective source S^g is generally divergent far from the worldline. To mitigate this unwanted 
behaviour, Vega et al. [41] multiply the 3-l-lD puncture field by a smooth windowing function, which 
has the effect of attenuating the effective source far from the worldline. Lousto and Nakano [37j have 
designed a puncture field with an effective source which is well-behaved at infinity, but is rather 
complicated to compute. In this work, we prefer to make use of a sharply-defined 'worldtube', 
which we introduce in the 2-|-lD domain after decomposition in azimuthal modes. 

The key idea was described in [59]. A worldtube T with boundary dT is constructed in the 
2-l-lD domain, to enclose the worldline (note that a 2-|-lD tube may also be interpreted as a 3-l-lD 
which spans the full range of azimuthal angles). Outside the worldtube, we evolve the modes of the 
retarded field ^I^^ti governed by the homogeneous wave equation. Inside the worldtube, we evolve 
the modes of the residual field ^>^, governed by the inhomogeneous wave equation sourced by the 
effective source modes 5"^^™ [defined as in Eq. (g]. Across the boundary of the worldtube, one may 
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FIG. 1: Illustration of a typical effective source derived via Eq. (|43|) from the 4th-order puncture field 



[given in Eq. (63)]. The plots show M^q ^^'g, for a particle on a circular orbit at ro = 7M on 



Schwarzschild spacctime, as a function of the coordinate differences Sr, 69, Sip, with the worldline being 
at 6r = 69 = S(p = 0. The source is shown on the spatial slices (clockwise from top left) 6(p — 0, 69 — 
and 6r ~ 0, with t = const. As expected, the 4th-order source is continuous but not differentiable at the 
worldline. 



convert between and $^ using the m modes of the puncture field, i.e. = + ^Ip . 

To summarise, 

= 5"" - n™^^ = S^, inside T, 

= 0, outside T, (44) 

= $™ - across dT. 

Here D'" is the d'AIembertian in the 2+ ID domain, obtained by making the replacement d^'/dcj)^ — >■ 
(—im)^ in the 3 + ID d'AIembertian □. 

For circular orbits, it is simplest to construct a worldtube with fixed coordinate widths F,. and 
Tg in the r and 9 directions. A worldtube of this form is illustrated in Fig. [2} 



III. IMPLEMENTATION: CIRCULAR ORBITS IN SCHWARZSCHILD SPACETIME 

Having laid out the general principles of the "puncture method with m-mode regularization" in 
Sec. [HJ let us move on to describe a simple implementation for the special case of circular geodesic 
orbits in the Schwarzschild spacetime. 
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FIG. 2: Visualization of the worldtube in 2+lD domain for a particle on a circular orbit. Here the t axis 
runs vertically, and the worldtube T is shown as a boxed domain of fixed width {Tr,Te}, centred on the 
worldline 7 at fixed r — rQ, 9 — tt/2. Inside the tube T we evolve outside the tube we evolve $™ti 
across the tube boundary dT we convert between the two using $™j = $^ + <f>p . 



A. Physical setup 

We consider the case of a pointlike test particle endowed with scalar charge q moving along a 
circular geodesic orbit at radius r = rg on a Schwarzschild black hole background. We work in the 
standard Schwarzschild coordinate system {t,r,6,ip}. Let z^{t) denote the particle's worldline, 
parameterized by proper time r, and let = dz^/dr denote the tangent vector. Without loss of 
generality, we assume that the orbit is in the equatorial plane, so that 

^^(.r) = [tp{T),ro,7T/2,ujtp{T)], (45) 

u'' = -^(1,0,0,0;), (46) 
Jo 

where tp^r) is the coordinate time on the worldline, u = (M/rj])^/^ is the angular frequency with 
respect to coordinate time t, and £ is the specific energy given by 

£^-ut = fo (1 - 3M/ro)-^/' , (47) 

where / = 1 — 2M/r and /q = f{ro). The retarded field $ret is governed by Eq. Q, explicitly, 

a^r,t = {-9r'^'\{-g)'^'g'"''^f] =S{x), (48) 

where g'^'^ = diag[— /"^, /, l/(r^ sin^ ^)] is the contravariant Schwarzschild metric and g = 
— r'^sin^^ is the metric determinant. The source term is defined by Eq. ([2]), explicitly, 

S = -^J^ir - ro)S - ^) 5 (^ - u:tp) . (49) 



B. Puncture scheme 



To briefly recap Sec. IIC-IIE a puncture scheme involves the introduction of a puncture field 

\n] fnl 

^-p , given analytically, whose singular structure is similar to that of $ret- The residual field 
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72 \TI\ 

is found from the difference between the fuU field and the puncture, i.e., = <I>ret — ^-p , and the 
effective source Ses is found from the d'Alembertian of the puncture via Eq. (43). After m mode 



decomposition, a worldtube (Fig. [2]) is constructed around the worldhne whose dimensions {Fj., Fg} 
are kept as controhable numerical parameters, of order M. The field equations to be evolved are 



given in (44) 



To proceed, we now require expressions for the m-mode decompositions of the puncture field 
and effective source, i.e., and S^. In this section we give explicit expressions for in the 
Schwarzschild spacetime for the 2nd, 3rd and 4th-order (n = 2, 3, 4) schemes, and we describe how 
to obtain and S^. 

1. Ist-order puncture, n — \ . 



A first-order puncture scheme was described in |59j . The first-order puncture field is simply 



= llHn ' ^here e[i] = [PrM^ + PeeSO^ + P^^5^'') ' , 



(50) 



with coordinate differences br = r — r^^bQ = Q — 7r/2, and Jc/? = 99 — as defined in (27). Here 
the coefficients are 



p 

± 7>' 



/o 



-1 



'^0' ~ ''"0 



2 ro - 2M 
ro - 3M ■ 



(51) 



Unfortunately, as discussed in Sec. II F, a first-order scheme is not sufficient to extract the SF; let 
us therefore proceed immediately to consider higher-order schemes. 



2. 2nd- order puncture, n = 2. 

A 2nd-order puncture scheme was described in [21 [60] . A 2nd-order puncture is given by 



[2] 
V 



Q/^[2], 



where 



with components 



,2 
£[2] 



eh + br {Qrrbr^ + QeeSO'^ + Q.^^bip'^) 



Qr 



M 
^0/0 



Qee = ro, Q 



ro - M 
ro - 3M 



(52) 
(53) 

(54) 



Within the m mode scheme we are motivated to re-express our puncture in terms of analytic 
periodic functions of b^p. Smoothness across = — vr, vr can be achieved by making the replacement 



bLp^ ^2(1- cos b^) = <5y?2 + 0(5(^4 



(55) 



in Eq. (53). This replacement does not affect the order of the 2nd-order puncture. Note that the 
alternative replacement bip^ — t- sin^ b^p is unsuitable, because it leads to an additional zero in e[2] 
at b^p = IT. 



The m-mode decomposition of the puncture field (52) and the effective source [obtained via 



(43)] is described in Appendix [B] With the replacement (55), it turns out that the m- modes can 
be found explicitly in closed forms involving elliptic integrals. 
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3. 3rd- order puncture, n — Z. 



Expressions for the 3rd and 4th-order punctures can be obtained using the covariant expansion 
method of Wardell and collaborators [651 ESI EO] • A 4th-order puncture is also employed in [H] . 
A 3rd-order puncture field is given by [HT] 

<l>?'=i^ + ^V (56) 



^[3] eW^fs], 



where 



M [dr^ + rlfo{59^ + 5^^)] [(2ro - 3M)r^'6r^ - rlfo{59^ + M^^] 
"^^1 - 6roVo^(ro-3M) ' ^''^ 

and 

3 3 

4] = 4] + E E U.,{S^if{6xjf. (58) 

i=l j>i 

Here the indices i, j run from 1 to 3, and we use the shorthand 6xi = 6r, 6x2 = 59 and 6x3 = 5(p. 
The coefficients are 

_ (8ro-M)M __r|^ rg/o(ro + M) 

" 12r4/3 ' " 12 ' ^^'^ " 12(ro - 3M) ' 

_ M(5ro - IIM) ro(3ro-2M)(ro-M) 
6ro/o' ''^^-6ro/o(ro-3M)' ""^^ " 6(ro - 3M) ' ^^^^ 

Again, we may replace the azimuthal variable to obtain a periodic function that is smooth 
across 5(p = —tt, tt. We make the replacement 

5^^ ^^-^ cos Sf + l; cos 26^ = 6^^ + 0(5/), (61) 
which preserves the order of the puncture. 



Next we compute the m-mode decomposition of the puncture field via (32), to obtain 

$^'™((5r, 69) = / $1^' {6r, 69, 6^') cos{m6^')d{6<^'). (62) 

27r 

Note that here we have used ip = 6(p + uot to factor out the time-dependence, and we have also 
used the symmetry of the puncture under 6ip ^ —6ip to eliminate the imaginary (sine) component 
of e~'™^'^, leaving the real (cosine) component. Analytic integrals are not readily available for 



computing (62) and so we resort to numerical methods. Thanks to the simple factorization of the 



t dependence in (62 ), we do not need to perform new integrals for each value of t in the simulation; 
nevertheless, we do need to compute the integrals numerically for every value of 6r and 69 within 
the worldtube. In future, it may be possible to find a way to obtain analytic representations for 
the m modes, by casting the puncture into a more tractable form. It is likely this question will 
require further investigation when eccentric orbits are considered (in which case the factorization 
of the time dependence is not so straightforward). 



The effective source S^g is found from inserting the d'Alembertian of (56) into Eq. (43). We 
used a symbolic algebra package to help with this calculation. The m-modes 5^ are obtained via 
numerical integration, in a similar manner to the above. 
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^. 4^th-order puncture, n = 4 



Finally, a 4th-order puncture field is given by [6T 



1 "[4] , ^[4] 



3 ~l~ 3 
^[4] £[2] 



(63) 



where 



3 3 

ef4] = 6f3] + M-i^r ^ ^ (5x,)2(5, 

i=l jr>j 



and 

Vrr 

Vre 



(6r§ - 2Mro + M2)m2 



(ro - M)M 



12 



(r§ + 4Mro - 9]VP)M 



_ (rg - 5Mro + 8M'^)M^ 
12rM' 12r3/2(ro-3M) ' = 

The remaining quantities are 

3 3 



12(ro - 3M) 
(3rg + 2Mro - 3M^)M 
6(ro - 3M) ■ 



[4] = "[3] + 7r:2 



6r2(ro -3M ^ 



with 



2(2ro - 3M) 



X„ = _ "v-u^ --/ ^ ^ -M-Vg(5ro - 3M), = -Pr\5ro - 9M)rl 



(64) 

(65) 
(66) 

(67) 
(68) 



Xre — 
and 



(2rg - 3Mro - 3M^ 
Mro/2 



2rl - 7Mro + 7M^ 



= -2M-V^(5ro - 6M),(69) 



3 3 



= 8ro2(!!-3M)^^'^^^-^'^"^^'^^"^' 



with 



Y 



2ro - 3M 

(rp + M) 
M/o ' 



= M-V^(3ro - 2M), Y^^ 
{rl - 12Mro + 21M2) 



M{ro - 3M)/o 



r|]/o(3rg -24Mro + 41M2) 
M(ro - 3M)2 ' 
2r^{3rl - 16Mro + ISM^; 
~ M(ro - 3M) 



(70) 

(71) 
(72) 



As in Sec. Ill B 3 , we used the replacement ( 61 ) and found the m-mode decomposition by performing 
the relevant integrals numerically. 



C. Simulation details 

In the following sections we employ the 4th-order puncture scheme unless otherwise stated. 
Therefore we will generally suppress the 'order' index [n] in the following sections, unless required 
for disambiguation. 
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1. Modal equations 



As discussed in Sec. II H the individual m- modes of the field are governed by a set of equations 
given in Eq. ( 44 ) . To make use of Eq. ( 44 ) we first require an explicit expression for the 2 + ID 
operator D'" on the Schwarzschild spacetime, i.e. 

□m^m ^ + f^rn^ + 2r-^{r - M)$^ + r^^ + cot 6*$^) - sin-2 6/$™. (73) 

Here can represent either or 'I'^ti depending on whether we are inside or outside the 
worldtube. Let us note that stationary retarded solutions of Eq. (48) fall off as 1/r towards spatial 
infinity, i.e. ~ 1/r as r — )• oo. This motivates the introduction of new field variables, 



The modes and are governed by the set of equations 



(74) 



f D^*^ = -(/r/4)5,^, inside T, 
0, outside T, 

across 9T, 



'—'it ^ ret 

— ^ret ' 



(75) 



where we recall that T and dT represent the interior and surface of a worldtube in the 2+lD 
domain, illustrated in Fig. [2j Here, 

= ^rn^_l_ [^m + _ ^2M/r + csc^ 0) , (76) 

where n and v are retarded and advanced Eddington-Finkelstein null coordinates, given by 



with the tortoise coordinate 



t — 



r + 2M In 



- 2M 
2M 



(77) 



(78) 



2. 2+lD grid and worldtube construction 

We construct a grid over coordinates u, v, 9, with linear spacing h in the u and v directions, 
and linear spacing A in the 9 direction. See Fig. [3] for an illustration. The grid may be seen as 
a stack (in 9) of causal diamonds (in u and v). The central r* = const line of each diamond is 



taken at the orbital radius, t^kq = '^*('^o) via (78). The two initial surfaces, at u = Ui = — r*o and 



V = Vi = r=KO intersect at t = 0, r = tq. The final surfaces a,t u = Uf and v = Vf intersect at 
t = ^maxi f = '''o- The shape of the grid, and use of null coordinates u and v, eliminates the need 
for boundary conditions except at the poles = 0,7r. 

Now let us introduce a worldtube of fixed coordinate widths {F^, , Fg}, centred on the worldline 
at r* = r*o, 9 = it/2. Consider an arbitrary grid point with coordinates (r*, 6). If |r=i, — r=i,o| < Fr, /2 
and \9 — tt/2\ < Tg/2, then the point lies within the worldtube; otherwise it lies outside. For 
convenience, we choose the worldtube widths F^^ and Tg to be integer multiples of the grid spacings 
h and A, respectively. 
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hy~h 




FIG. 3: 2+lD finite difference scheme. The left plot shows the grid in u, v and 9. The right plot shows a 
single 'cell' for the finite difference method. These plots are reproduced from [SI], Figs. 1 and 2. 

3. Initial and boundary conditions 

For each m, we specify a 'zero' initial condition on the initial surfaces u = Ui and v = Vi, 

^Tedui, V, e) = 0, ^T,,{u, v„ e) = 0. (79) 

The initial condition is not a solution of the sourced field equations. However, we expect that 
'junk' in the initial condition will radiate away, towards the horizon and infinity, and that after 
sufficient time the field near the particle will approach the correct (retarded) stationary solution. 

At the poles, boundary conditions are required. As argued in ||59j, the physical boundary 
conditions to be applied at the poles are 

de^T,T\0 = Cvr) = 0, M/^f^e = 0,vr) = 0. (80) 

To implement these conditions, we simply set = at the poles for m 7^ 0, and for m = we 
extrapolate to obtain 

*™=°(^ = 0) = 3 [4*"r°(^ = A) - "^TefiO = 2A)] + O (A^) , (81) 
^^='{0 = = I [^^TefiO = vr - A) - ^TcfiO = ^ - 2A)] + O (A^) . (82) 

Here the error term is 0{A^), rather than 0{A^), since ^'^^ is an even function of ^ at = (and 
an even function of vr — ^ at = vr). 

4-. Finite difference method 

We employed the finite difference method which was described, applied and tested in 
Figure [3] shows a 'cell' of grid points with centre c. Let us assume (for now) that all the grid points 
in the cell lie outside the worldtube so that ^'j", . . . , represent values of the retarded field mode, 
at the positions shown. Finite-difference approximations for the values of the field and its 



19 



derivatives at c are obtained through the following expressions: 



/l2 



2 (^-^ + ^'™) 



2A2 



+ (A^ 



c 



4A 

^2 +^3 



+ O (A^) 



(83) 
(84) 
(85) 
(86) 



Now, assuming that the values at points 2 to 8 have been obtained in previous steps, we may insert 



(83)-(86) into Eq. (75) and rearrange to find 



^ [("f ^ + + ^-7* + - 2*^ - 2*^) /A^ + cot 9 + ^^-^^- ^^) /(2A) 
- (2M/r + csc^ 6) (^^ + ^'^)] + O {h^A\ /i^) . (87) 

Here, all r, 6 dependent coefficients are evaluated at the centre of the cell. Conversely, if all points 
in the cell lie inside the worldtube, then . . . , ^™ represent values of the residual field mode, 
^jl, and one may repeat the argument above to obtain 

(88) 



= [RHS of Eq. (87)] + h^Z, 



,2 rym 



where Z, 



cff 



-frSTs/4. 



In the case where the finite difference cell straddles the boundary of the tube (so that some 
points are 'in' and some 'out'), we make use of the puncture field in the following way. If point 
1 (Fig. [3]) is 'out', then we first demote all 'in' points in the cell to 'out' points using <I>™ = + 
before applying (87). If, conversely, point 1 is 'in' then we promote all 'out' points in the cell using 



$P before applying (88). This strategy is discussed in more detail in Sec. VB of 
For a fixed ratio A/h, the finite difference equation (87) has a local discretization error of 
0{h'^). In vacuum, therefore, we expect (and find) the scheme to be quadratically convergent 
[i.e. to exhibit a global accumulated finite-differencing error which scales as O^h?)]. Unfortunately, 
quadratic convergence of our simple scheme is by no means assured if a non-smooth source term 



is present, as in Eq. (|88|). Consi der th e special case of a grid cell whose centre c lies exactly on the 

II h[ for puncture orders n < 4 the effective source 5'^^^'^' is not 



worldline. As discussed in Sec. 



continuous across the worldline, and Sg^^^'™ cannot easily be evaluated for this cell. A strategy 
for dealing with this problem was discussed and implemented in |59] . The case of a cell centred on 



the worldline was handled separately, taking into account the singular structure of S^^^ . A similar 
method was employed here, in the cases n = 2 and 3. The procedure leads to a local error in the 
central cell that scales as h? In /i, which leads to a term in the global accumulated finite-differencing 
error that scales with In h. Although this undesirable behaviour could perhaps be eliminated 
by using a more sophisticated finite-difference scheme in the vicinity of the worldline, we have not 
pursued such an approach here, principally because we now have a 4th-order puncture {n = 4) 
available. In the 4th-order case, no difficulty is encountered for the cell on the worldline (since 
the source S'^q* is continuous across the worldline) and the global convergence rate is found to be 
quadratic. 



5. Numerical stability 

Finite difference methods can suffer from numerical instabilities, where generic small-amplitude 
short-wavelength perturbations (originating for example from truncation errors) are amplified ex- 
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ponentially, eventually overwhelming the physical solution. Our finite difference method suffers 
from a numerical instability if the ratio of grid spacings A//i is set to be too small. We observed 
in our vacuum simulations that this instability arises first near the poles, and appears as a spu- 
rious oscillation in the direction with wavelength 2A and an exponentially-growing amplitude. 
In Appendix [C] we apply a von Neumann stability analysis (see, e.g., [71j) to our finite difference 
equation (87) which suggests that a necessary condition for stability is 

^ > ^max (r- 7^/2^ ^Jl + m?/A , (89) 



where max(r ^/^^^) ~ 0.19245M ^. Equation (89) becomes a highly restrictive condition when 
m is large. For m > 3 this condition becomes stronger than the standard Courant condition, 

^ > max (r- V^/^) , (90) 

which is obtained by insisting that the numerical domain of dependence contains the physical, 
continuum domain of dependence at each point in the evolution. 

To mitigate the instability in large-m modes we may move the grid boundary inwards from the 
poles. First we note that solutions have a simple asymptotic form near the poles, i.e. 

"^TedO < 1) = AO"^ + ^^"+2 + O(0™+^), (91) 
*JSt(^ -0-^1) = A{tt- er + B{7r - 0)^+2 + 0{{7r- 0)""+^) , (92) 

where A and B are constants. In other words, the large-m modes are very 'flat' near the poles. If 
we move the boundary point inwards from 9 = to 9 = kA (with k being a small positive integer) 
then, repeating the analysis of Appendix O the stability condition becomes 



A 1 

— > - max 
h - 2 



(r-i/'/') Vl + "iV[4(A: + l)2], (93) 



which is less restrictive than the original condition (89). The boundary condition near 9 = 
changes to 

k = l : ^"^^^^9 = A) = - [2=^-™^™ (2A) - 3^-'"^St(3A)] + C'(A"^+^), (94) 

k = 2 : ^'"t^^e = 2A) = ^ [12 X (2/3)™ ^™t(3A) - 5 x 2"'"^'™ (4A)] + 0(A"^+^), (95) 

etc. The modes are symmetric under 9 ^ it — 9, so equivalent boundary conditions may be applied 
near the south pole at 9 = tt. In our implementation, we fix the ratio A/h and, for a given mode 
m, we find the minimum value of k required for stability using (93). 



D. Simulations and data extraction 



with boundary conditions (80) 
dimensions Uf — t 



For a given m, we evolve the initial data (79) according to the finite difference scheme (87) 



Vf 



Vi 



([82]) [or (l94|-(|95j)] on the 'diamond stack' 2+lD grid (Fig. |3[ of 
^maxi to obtain a numerical estimate for an m-mode residual field, 
^'^(t, r, 6). Of course, the numerical solution obtained depends on the 'physical' parameters, tq/M 
and m, and a set of 'numerical' parameters, {num.} = {/i. A, F^*, Fg, tmax, • • •}• We will refer to a 
simulation for a particular m, rg/M with a unique set of {num.} as a 'run'. 

As described in Sec. II G , the SF and the R field (at a worldline point z) are computed from a 



sum over modal contributions -F™ and $7^(2) [see Eq. (35) and (34)]. Let us briefly describe how 
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the modal contributions are obtained from our 'runs'. First we assume that the total time tmax is 
sufficiently large that the residual field inside the worldtube has settled into a quasi-stationary state 



at late times (see Sec. IV A 5| for further consideration). Next, we read off the following quantities, 



n 



(ti,ro,7r/2,wti) , 



n 



(ti,ro,7r/2,a;ti), 



2mroilm[^'?? (ti,ro,V2) e 



(96) 
(97) 
(98) 



where 



ro$7e with as defined in Eq. (33). Here all quantities are evaluated at grid points 
on the worldline at a late time t = ti, where ti < tmax, and the derivative with respect to r* is 
found via central differencing on the grid. The SF and radiative field on the worldline are found 



by inserting (96)-(98) into the mode sum reconstruction formulae (35) and (34), respectively. In 



practice, we compute ^r, F^'^^^ and F^^^^ directly from mode sums; the remaining components are 



pself 



given by F^ 



'self 



-UJ 



F^^^^ and, by symmetry, = 0. 



E. Sources of numerical error 



Of course, the values extracted from a particular run via Eqs. (96)-(98) depend in part on the 



set of numerical parameters {num.}. Inevitably, the values contain numerical error, which we may 
define as the difference between a given numerical solution and the (unknown) exact solution. To 
compute accurate SF estimates we must first seek to understand the various sources of numerical 
error that arise in our implementation. By judiciously combining the results of multiple runs, we 
then attempt to quantify and minimize the error. 

Let us identify several key sources of numerical error. These will be more fully described in the 
next section. The evolution of a single mode is affected by the following: 



Discretization error (Sec. IV A 2), associated with use of a finite grid spacing /i. A, i.e. 



0, A ^ 0}). 



(99) 



Worldtube error (Sec. IV A 3). Changing the dimensions of the worldline {rriTg} affects the 
amplitude of the discretization error, but should not affect its scaling with h. 



Source cancellation error (Sec. IV A 4), associated with roundoff error arising in the calcu- 



lation of S'eff close to the worldline, from the delicate mutual cancellations of large terms in 
the high-order puncture. 



Relaxation time error (Sec. IV A 5), associated with the time it takes for junk radiation to 
decay, and the solution to reach a steady state, i.e. 



oo 



(100) 



In computing mode sums, there arises further errors: 



m-mode summation error (Sec. IV B 2). Only a finite number of modes may be calculated 
numerically. We impose a large-m cutoff /n-max, and estimate the contribution from the 
remaining modes by fitting an appropriate model. 



Mode cancellation error (Sec. IV B 3). If the magnitude of individual modal contributions 
(or $?^) is large in comparison to magnitude of the total mode sum F^^^ (or ^>r), then 
the relative error in the mode sum may be much larger than the relative error in individual 
modes. 
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r*/M 9 

FIG. 4: Field modes on constant time slices (at t = tmax/2) for a circular orbit at vq = 7M (r*o ~ 8.8326M) 
with the 2nd-order puncture scheme. The left plots show field modes at fixed 9 — Tr/2 and the right plots 
show field modes at fixed r — vq, for a range of modes m = 0, 1, 2, and 5. Inside the worldtube (visible as 
the central 'trough'), the dashed (red) line shows the full retarded field and the solid (blue) line shows 
the residual field, 'I'^. The numerical parameters are {h = M/8, A = 7r/40,r9 = TT/2,Tr, = 5M}. 

IV. RESULTS AND ANALYSIS 

In this section we present a selection of results from our numerical simulations. We discuss the 
challenge of minimizing numerical errors by giving illustrative examples. 

A. Individual m-modes 

1. Simulations and visualisation 

Let us consider a typical 'run', i.e. a simulation for a single m and rg and a unique set of 
numerical parameters {num.}. The results of a run can be visualised by examining particular 
slices through the u-v-9 grid. Three informative slicings are: (i) t = tmax/2, = vr/2, i.e., across 
the central line of the uv diamond in the equatorial plane, (ii) t = tmax/2, = r^o, i.e., from pole 
to pole, and (iii) r* = r*o, = 7r/2, i.e., 'along the worldline'. 

Figure [4] shows m-mode contributions to the field modes along the constant-i slices (i) and (ii), 
for the 2nd-order puncture scheme. The worldtube is visible as the 'trough' in the centre of these 
plots. The residual field (solid line) is continuous and differentiable across the worldline (at 
6 = 7r/2, r = ro), whereas the retarded field (dotted line), found using ^'^^ = + ro$p, 
diverges at the worldline. Figure [5] shows m-mode contributions to the field modes 'along the 
worldline' as a function of time, i.e., on slice (iii). At early times, the signal is dominated by 'junk 
radiation' arising from our imperfect choice of initial condition. The effect of the junk radiation 
diminishes with time, and the field approaches a steady state. 

Data on constant-f surfaces can alternatively be visualised using 3D plots. Figure [6] illustrates 
field modes as functions of r^, and 9, on constant-t slices. Here, the worldtube is apparent as a 
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FIG. 5: Modes of the residual field, ^^tj, as a function of time, evaluated along the worldline = r.^,o, 
9 — tt/2 [slice (iii) in text]. The initial burst of 'junk radiation' (due to the imperfect initial condition) 
radiates away, and the field approaches a steady-state in the vicinity of the worldline. 



thin interior rectangle. The central feature becomes sharper as m is increased, with the solution 
becoming 'flatter' at the poles (i.e. near ^ = and vr). The oscillations seen in Fig.|6]at large r are 
outgoing waves emitted by the particle; they have a wavelength of ~ 27r/(ma;). 



2. Discretization error, convergence tests and Richardson extrapolation 



Figure [5] illustrates how the system approaches a steady state once initial 'junk' radiates away 
(we return to consider this point more carefully in Sec. IV A 5). In the steady-state regime, we 



may extract estimates for the m-mode contributions on the worldline, Eqs. (96)-(98). The values 
obtained obviously depend upon the grid spacings, h and A. We set /i to be a simple fraction of 
M, i.e. h = M/rij-Qs with n^es an integer. Then, rather than varying {/i, A} separately, we fix the 



ratio A/h. Here, we are limited by the stability condition (93) which imposes an (ro-independent) 
constraint upon A/h. For convenience we choose A, h such that Ores = /ivr/(MA) is a fixed integer 
(typically arcs = 10), and we determine k (i.e. the displacement in grid points of the numerical 



boundary from the poles) for each m according to stability condition (93). 

The left plot of Fig. [7] illustrates the (4th-order) residual field as a function of time, for a range 
of resolutions rij-es = 16, 24, . . . 64. It suggests that the field converges towards a limiting curve as 
-)• oo. We may test the convergence rate by taking ratios of the results of runs at different 



n 



resolutions n^es- For example, consider the ratio 



Xih) 



X{4h)-X{2h) 
X{2h)-X{h) ' 



(101) 



where X G {^tI^ F^} and X{kh) denotes the extracted result from a run with grid spacing kh. 
If the convergence rate is quadratic (i.e. if the dominant term in the numerical error scales as /i^), 
then X would approach the value of 4 as /i —t- 0; on the other hand, if the convergence is only linear 
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FIG. 6: Sample numerical results for the 4th-order puncture scheme at rg — 7M. Here we show (outside 
the worldtube) and 4'^ (inside the worldtube) as a function of and 9 at late time t = 200Af for m = 0, 
1, 5 and 10. The worldtube is visible as a thin rectangle of fixed width Tg, Tr* around the particle position 
at = 7r/2, = r*(ro) ~ 8.83M. 



then we expect x ~^ 2. In Table |ll| we present sample values of x for both ^IJ^ and for the 2nd, 
3rd and 4th-order puncture schemes. The data in Table |ll] shows convincingly that the 4th-order 
scheme is quadratically convergent (i.e. x ~^ 4). The data for the 2nd and 3rd-order punctures 
is less conclusive. In {59j it was noted (for the 1st order puncture scheme) that the irregularity of 
Scs at the worldline disrupts the global quadratic convergence of our finite difference scheme. In 



Sec. IIIC4 we described how the procedure for evaluating in cells on the worldline is expected 
to introduce an additional term in the global discretization error which scales with h^lnh. To test 
for the presence of this term, we construct the ratio 

X{8h)-5Xi4h) + 4X{2h) 
XiogW - ^^^^^ _ ^^^2/,) + 4^x{h) ' > 

with xiog — ^ 4 as /i — >• 0, if a /i^ In term is present. Table [ll| gives convincing numerical evidence 
in favour of the presence of an /i^ In/i term, at 2nd (field and SF) and 3rd order (SF only). 

In order to improve our estimates of the 'physical' results, we used a fit model to extrapolate to 
/i — >• ("Richardson's deferred approach to the limit" [H]). As discussed above, the appropriate 
fit model depends on the order of the puncture. We use Xq + Ah'^ + Bh^ Inh + 0{h^) for ^^"^^"^ 
and i7;!"=2'3l'^^ and Xq + Ah^ + 0{h^) for ah other cases. The fitting procedure is illustrated in the 
right panel of Fig. The example shows that the 4th-order data for a field mode is well fitted by 
the simple model ^'?g(^) = ^'™(/i = 0) + Ah"^ + Bh^ . 
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FIG. 7: Test of convergence. The left plot shows a mode of the residual field on the worldline, ^'7^, as a 
function of time, for resolutions Ures = 16,24,32,48,56, and 64 [here h = M/rires, A = 7r/(10ni.es), m = 2, 
ro = 6M, with a 4th-order puncture]. The right plot shows the field extracted at late times (t = 300M) 
as a function of grid resolution. The line shows the best fit extrapolation model, (—1.07487 + 1.08598/i^ — 
1.0056/1^) X 10-2. See also Fig.^ 
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X 


Xlog 


X 


Xlog 


X 


X 


m — 
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m = 
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3.30 
3.24 
3.13 


4.05 
4.16 
4.23 


3.92 
3.95 
3.90 


6.62 
8.87 
7.15 


3.85 
3.97 
3.87 


3.98 
4.00 
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7M 
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Xlog 
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Xlog 
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X 


m ~ 
m — 
m = 


5 

10 
15 


3.39 
3.08 
2.90 


4.07 
4.11 
4.16 


4.25 
4.72 
0.82 


3.86 
3.87 
4.10 


3.96 
3.92 
3.86 


3.99 
3.99 
3.98 



TABLE II: Sample convergence tests for ^f^' and i^™ for the 2nd, 3rd and 4th-order puncture schemes. The 
ratios x xiog are defined in Eq. (101) and (102). Here x is an additional convergence ratio using points 
closer to the asymptotic regime, defined by x = (135/52)[F;"(nros = 48) - F^'in.^s = 56)] / [F^ (ures = 
56) — i^™ (fires = 64)]. For the 4th-order puncture, the ratios x and x are ^ 4, implying that convergence is 
quadratic. Similar behaviour is apparent for the residual field (but not the SF) at 3rd order. In all other 
cases shown {'^n and F™ at 2nd order and F™ at 3rd order) quadratic convergence is not clear. The data 
suggests that xiog ~ 4 , which implies that global convergence is affected by the presence of an 0{h^lnh) 
term, due to the non-smoothness of the effective source on the worldline. Similar results are found for the 
angular component F™. 



3. Worldtube error 

It is important to check that the dependence of the numerical results upon the dimensions of 
the worldtube (i.e., F^* and Tg) diminishes as rires — ^ 0. Figure [s] shows results from using three 
different worldtubes: narrow (F^* = 1.25M, Vq = vr/8), medium (F^* = 2.5M, Tq = vr/4) and wide 
{Tr* = 5M, F51 = 7r/2), as a function of grid resolution. The plots show that the magnitude of the 
discretization error (but not its scaling with h) depends on the worldtube dimensions. We would 
expect the magnitude of the discretization error to scale with the maximum absolute value of the 
numerical variable, which is typically the value of the field mode just outside the worldtube (see e.g. 



26 



-1 .04e-02 
-1.05e-02 
-1.05e-02 
-1.06e-02 
-1.06e-02 
-1.07e-02 
-1.07e-02 
-1.08e-02 
-1.08e-02 



M ' 

N3rrow 


1 1 


1 1 1 


Medium 






- Wide 


□ 








/ 












/ 

/ 






/ 


- 




y 

/ 






/ 




























^^^^ 








, - -X - 






- -X - - " 















6 ^ 



0.01 0.02 0.03 0.04 0.05 0.06 

h / M = 



-4.48e-04 i r 

Narrow 

-4.49e-04 -Medium 
Wide 

-4.506-04 - 

-4.516-04 - 

-4.526-04 - 

-4.536-04 - 

-4.546-04 - 

-4.556-04 - 

-4.566-04 - 
-4.576-04 
-4.586-04 







0.01 



,02 0, 
h/ 



0.04 0.05 



0.06 



FIG. 8: Finite worldtube size effect and extrapolation to zero grid spacing. Tlie plots show typical (m = 2, 
To — 6M, 4th-order puncture) modal contributions to the residual field (left) and radial SF (right) as a 
function of grid resolution (h = Af/rtres, c^res = 10), for worldtubcs of three different widths, i.e. (i) narrow: 
T„ = 1.25M, Te = tt/8, (ii) medium: F^* = 2.5M, Tg = 7r/4, and (iii) wide: F„ = 5M, Tg = 7r/2. 
Numerical results from runs at six resolutions are shown {h = M/n,-cs where n^^s — 24, 32, 48, 56, 64 
and ttros = 10) as data points. The best fits to the model Xq + Ah? + Bh^ are shown as lines. The 
extrapolated values for the field (radial SF) vary only within ~ 0.0004% ['-^ 0.0005%), which is considerably 
less than the relative error of ~ 0.024% (0.036%) obtained by comparing the highest resolution result with 
the extrapolated value. 

Fig.|4]). Hence we expect that using a wider worldtube will generally decrease the magnitude of the 
grid resolution error, and this is what is observed in Fig. [8j Note that arbitrarily large worldtubes 
are not practical however because (i) the computational expense of calculating S*^ scales with the 
spatial cross-section of the tube, and (ii) generally diverges far from the worldline, diminishing 
numerical accuracy. 



4-. Source cancellation error 

The 4th-order effective source is continuous across the worldline, and is zero on the worldline. 
It is calculated from the d'Alembertian of the puncture field, which is divergent at the worldline. 
The calculation of Seff near the worldline involves the delicate cancellation of large terms. This 
calculation is susceptible to numerical round-off error. 

We found that, unmitigated, source cancellation error has greatest relative impact on high- 
resolution runs (which have a greater density of grid points in the vicinity of the worldline) , and 
the error disrupts the smooth convergence to infinite resolution exhibited in Fig. [7] and Fig. [8} 
In consequence, unmitigated source cancellation error affects the validity of the extrapolation 



described in Sec. IV A 2 To deal with the problem, we used a symbolic algebra package (Maple) 



to obtain an approximation for the source close to the worldline. First we introduced a scaling 



parameter A (see Sec. II D) via = X6x^, and expanded the full expression for Sgg in powers of 
A at A = 0. We verified that (for the 4th-order puncture) the coefficients of the divergent terms 
(at orders A~^, A~^, A~^), as well as the constant term A*^, are identically zero, 

SesiXSxf") = Xsl{6x^') + ©(A^), (103) 

where si is a (i.e. discontinuous but bounded) function of rescaled coordinate differences 6x'^. 



Equation (103) is not sensitive to large numerical round-off errors in the vicinity of the worldline. 
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TABLE III: Sample data for relaxation error due to dissipation of junk radiation in low modes (to = 0, 1, 
2) for ro = 6A/ and ro = 20M. Here we give numerical data for A$™ = $™(t = 300M) - l>™(250M) and 
AF™ = F™(t = 300M) - $^(250M), i.e the difference between field and SF values extracted at t = 300M 
and t — 250M, after extrapolation to /i — has been performed. The data shows that (i) the magnitude of 
the error decreases with m, as expected from considering power-law decay; (ii) for to = 0, the relative error 
in the field mode, A^^^^/l)^"", is larger than the relative error in the radial SF mode, AF™="/F™"°; (iu) 
for TO = 0, the absolute error increases somewhat with radius rg; hence the relative error increases rapidly 
with rg (see text). 



Using Eq. (103) very close to the worldline, and the full expression further away, significantly 



reduces the effect of source cancellation error. 



5. Relaxation time error 



We wish to estimate the steady-state values for the mode-sum contributions. Obviously it is not 
possible to run the simulation for an infinite amount of time, and long runs are computationally 
expensive. With our 2-l-lD grid (Fig. |3|, doubling the physical simulation time (tmax — ^ 2tmax) 
quadrupoles the run-time (i.e. the CPU time). We explored a range of methods to obtain accurate 
estimates of steady-state values from simulations with finite tmax, which we briefly describe below. 
We begin with an estimate of the approximate magnitude of the errors. 

a. Magnitude of error. Table |III| provides data on the approximate magnitude of the relax- 
ation error in the lowest modes m = 0, 1, 2, for orbits at ro = 6M and ro = 20M. Let us define 

= $™(t = 300M) - ^^{250M), i.e. the difference between modal contributions 'read off' 
at t = 300M and at t = 250M, and define AFJP in a similar way. We make the following simple 
observations: (i) |A<I>^| and |AF™| decrease in magnitude as m increases, i.e. the dominant error 
is in the m = mode. This is expected, since the m = mode contains the monopole which 
relaxes most slowly (see below); (ii) for the m = mode, the relative error A^^^^/^^^'' is larger 
than the relative error AF^^^ /F^^^; (iii) the absolute errors A$J^ and AF^ increase somewhat 
in magnitude in going from ro = 6M to ro = 20M; consequently, the relative errors A<I>^/<I>^ and 
AF^/F^ are significantly worse at ro = 20M than at ro = 6M because the total field and SF 
diminish rapidly as rg increases. 

b. Power-law relaxation. The data in Table |III| suggests that the m = mode relaxes to 
equilibrium most slowly, as expected. Figure [9] shows the relaxation of the m = mode of the 
residual field (left) and the radial SF (right) as a function of time, for a range of radii. The field 
exhibits a power-law relaxation, i.e.. 



it) 



^^="{t ^ oo) + At-"^ + 0{t- 



(104) 



28 



>fS< 



lis, 



lis, 




400 600 800 1000 1200 1400 1600 1800 2000 



1.634316-03 - 




c 1.71686-04 



-=q 1.71676-04 



rg = 6M 



400 600 800 1000 1200 1400 1600 1800 2000 



rn = 1 0M 



400 600 800 1000 1200 1400 1600 1800 2000 




400 600 800 1000 1200 1400 1600 1800 2000 
t/M 



400 600 800 1000 1200 1400 1600 1800 2000 
t/M 



FIG. 9: Relaxation towards equilibrium of the m = modes for a range of orbital radii (rg = 6M, lOM and 
20M). The left plots show the m = mode of the field, 'i>^^^, evaluated on the worldline, and the right 
plots show the same mode of the radial SF, as a function of time t, for a low-resolution long-time 

run (uros = 16, a^cs = 10, imax = 2000 Af). In the appropriate late-time regime, the data are well-fitted by 



simple power-law relaxation models (see Fig. 10 ) 



It is expected from theory |72j that the appropriate index for the monopole component of the m = 
mode (i.e. for the / = multipole) is ry = 3 if the initial junk radiation is localized in space and 
r/ = 2 otherwise. In this case, the latter index is appHcable because the steady-state solution for 

oo. Figure 



tends to a non-zero value in the limit r 



10 



shows the numerically-determined 



local power-law index, defined by r/(t) = —t^j^^^/"^^^^ — 1, plotted as a function of time. For all 
orbital radii, the local index asymptotes to 2 in the late-time regime, as expected. Making use of 
this observation, we may minimise relaxation error in the m = mode of the field by fitting the 
numerical data to a power-law model (104), with rj = 2, to extract the steady-state value. Whilst 
this procedure is straightforward for the m = mode, it is more difficult for higher modes (m > 0) 
which also exhibit damped oscillations of frequency rrno. However, it suffices for our purpose to fit 
only the m = mode, since this is by far the dominant source of relaxation error. 

The relaxation of the m = mode of the radial SF also exhibits power-law decay, but in this 
case, the appropriate index is = 3 as shown in the right-hand plot of Fig. [Toj It turns out that 
the slowest-decaying part of the monopole (/ = 0) perturbation in ^^^^ does not depend on 
radius [see Ref. [73j, in particular Eq. (89)] and as a result the relaxation of F™^*^ is one power 



of 1/t faster than naively expected. The right-hand plot of Fig. 10 shows that the onset of the 
late-time regime, where power-law relaxation is manifest, increases with orbital radius. Unless one 
can evolve for very long time, fitting a simple power law to the radial SF generally does not give 
good results. Longer runs are computationally expensive, since the runtime and memory usage 
scales as t^ax- practice, the computational burden associated with high-resolution (rires ^ 64, 
Ores = 10), long-time (tmax ^ lOOOM) runs may be prohibitive. This leads us on to consider an 
alternative strategy. 



29 



i3< 



o 



O 
1-1 



2.3 



r' =6M --^ 

r„ = 7M 

rn = 8M 

ro= 10M 

rn = 14M - ■ - - 
rn = 20M -- -- 



II 



O 



o 



1200 
/M 



1400 1600 1800 2000 



4 
3.5 

3 
2.5 



rUeM 
To = 7M 
r(, = 8M 
ro = 10M 
rn= 14M 



1000 1200 
t/M 



1400 1600 1800 2000 



FIG. 10: Power law relaxation of the m — Q mode. The left-hand plot shows the local power-law index ij of 
the relaxation of the m = mode of the field, determined from ri{t) = — t^'^^''/\E'^=° — 1 (where overdot 
denotes differentiation with respect to t, and derivatives are evaluated numerically). The right-hand plot 
shows ri{t) for the radial SF mode the case of the field (left), the index tends towards = 2, as 

expected for a non-compact I = perturbation. In the case of the radial SF (right), the index tends towards 
77 = 3, although for large tq it takes a long time to reach this asymptotic regime. The higher index in the 
right plot (i.e. 77 = 3) is due to the fact that the slowest-decaying part of the monopole is spatially constant 




FIG. 11: Three-stage multigrid refinement method. A 'crude' run (1) provides initial data for an 'interme- 
diate' run (2) which in turn provides initial data for a 'fine' run (3). Here we show a single slice of the grid 
a.t 9 = const. 



c. Multigrid refinement. A complementary solution to fitting a power-law relaxation model 
is to use a kind of 'mesh refinement'. The key idea here is to run the bulk of the simulation at 
low resolution (which is computationally cheap) and then improve the resolution in the late-time 
regime. In other words, we use the results of low-resolution runs to improve the initial data used 
in the final high-resolution run, which in turn reduces the amplitude of the final relaxation error. 

Figure 11 gives an illustration of a three-stage process of grid refinement. Here the 'fine' grid 
(3) has twice the resolution (in both radial and angular directions) of the 'intermediate' grid (2), 
which in turn has twice the resolution of the 'crude' grid (1). A 'crude' run (1) provides a rough 
estimate for the field everywhere in the largest grid. Initial data for the 'intermediate' run (2) is 
then obtained by interpolating the values of the field read off along the initial boundary of grid 2. 
In a similar way, the intermediate run (2) then provides data for the fine run (3). The 'fine' grid 
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takes approximately twice as long to run as the intermediate grid, and four times as long as the 
crudest grid. It is much faster to run the multigrid scheme than to run the single large grid (1) at 
the highest resolution. The speed-up factor for a three-level grid is approximately 8^/(1 -|-2-|-4) ~ 9, 
and it is 8^/(1 2 4 -h 8) f« 34 for a four-level grid. 

Figure [12] shows typical results from a multigrid implementation, for orbits of radii tq = 6M 
and ro = 30M. Here, the results of multigrid refinement are compared against unigrid results 
(which, as argued above, take much longer to run). Junk radiation is visible in the 'intermediate' 
resolution which starts at t = 500M and in the 'fine' resolution which starts at t = 750M. After the 
high-frequency junk has dissipated, at late times, the multigrid and unigrid results are found to be 
in close agreement. Importantly, the difference between multigrid and unigrid results at the same 
resolution is observed to be much less than the difference between unigrid runs at different resolu- 
tions. In other words, these plots demonstrate that the left-over error associated with refinement 
(e.g. from the interpolation procedure) is much smaller than the discretization ('grid resolution') 
error, and that grid refinement is a useful technique which greatly diminishes the computational 
burden. The only constraint upon the scheme is that the finest grid must be large enough that the 
high-frequency junk arising from interpolation has time to dissipate. 

The simple method of multigrid refinement outlined here is crude in comparison with the more 
systematic adaptive mesh refinement (AMR) algorithm implemented by Thornburg [62] in 1-l-lD. 
Our hope is that an AMR scheme could be applied to 2-l-lD simulations in the future [74J. 



B. Mode sums 



Now let us turn attention to the mode sums given in Eqs. (34)-(35) and their numerical calcu- 
lation. 



1. Large-m asymptotics and convergence 



Let us consider the behaviour of the modal contributions 
limit 



rpm 

r 



and 



in the large-m 



As we argued in Sec. II G their limiting behaviour depends on the order of the puncture 
scheme. In this section, we give numerical evidence in support of the following conclusions: (i) 
~ 0(m~^) for 2nd-order punctures, and ~ 0{m~^) for 3rd and 4th-order punctures; (ii) 
~ 0{m^'^) for 2nd and 3rd-order punctures, and ~ 0{m~^) for 4th-order punctures; 
and (iii) F^,FJ^ oc exp(— /3m), where f3 is an m-independent positive constant which depends on 
orbital radius tq. A heuristic explanation for these behaviours was given in Sec. II G[ 

Let us examine the magnitude of the modal contributions for a circular orbit of radius ro = 7M, 



for implementations of the 2nd, 3rd and 4th-order puncture schemes. Figure 13 shows the modal 
contributions to the residual field, The upper plot shows that the modal contributions can 

change sign; in particular, the m = and m = 1 modes have opposite signs. The lower, log-log 
plot suggests the scaling ~ 0{m~'^) in the large-m limit, with C = 2 for the 2nd-order puncture 
and = 4 for the 3rd and 4th-order punctures, as anticipated in Sec. II G (see Table |l]). Figure 
14 shows modal contributions F^, for the conservative component of the SF. Here again we see in 
the upper plot that low-m modal contributions can take either sign. We also see strong evidence 
for power-law convergence, i.e. -F^ ~ 0{m~''), with C = 2 for 2nd and 3rd-order punctures, and 
= 4 for the 4th-order puncture. This, again, is consistent with the predictions of Sec. II G 
Figure 15 displays the modal contributions to F^^^^^, the dissipative component of the SF. The 
plot shows that the modal values are independent of the order of the puncture (up to numerical 
error), as foreseen in Sec. IIGl, Furthermore, the modes exhibit a clear exponential convergence. 
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FIG. 12: Sample results comparing multigrid and unigrid evolutions of the m — mode. The left-hand 
plots show the evolution of the field mode XP^^*^ [upper] and radial SF mode F™^" [lower] at rg — 6M. The 
right-hand plots show similar data at tq = 30Af . The dashed lines show unigrid evolutions at three different 
resolutions, n^cs = 4, 8 and 16. The solid lines show the results of the multigrid scheme. Refinements in 
resolution occur at t = 500M (rii-cs = 4 to rircs ~ 8) and t — 750M {rii-cs 8 to rtrcs = 16). After a period 
of transition dominated by junk radiation, we find that the 'refined' data asymptotes to the unigrid data. 
Note that the multigrid evolution is approximately 9 times faster than the unigrid evolution at equivalent 
resolution. 



oc exp(— /3m) (with /3 > 0) at large m. 

Finally let us examine the dependence of the modal contributions upon the orbital radius rg, 
focussing on the 4th-order puncture. Figure 16 shows the magnitude of various modal contributions 
to the radial SF for a range of radii. The magnitudes of the modal contributions diminish with 
increasing rg, though the relative contributions of different modes do not change substantially. The 
plot makes it clear that the 'asymptotic regime', in which the modal contributions follow an inverse 
power law in m, begins at around m ~ 10 for all radii. 
Figure 
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shows the magnitude of modal contributions to the angular component F^'^^^, 



for a 



range of radii. Again, the magnitudes of the modes diminish as tq is increased. Exponential decay 
of the modes with m is seen for all radii, with the decay rate f3 increasing with tq. 



2. Mode summation and large-m fitting 

For the field {^r) and conservative SF [F^^^^) it is important to account for the modes in the 
large-m 'tail', whereas for the dissipative SF {F^^^) this is not necessary (as the large-m modes 
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FIG. 13; Modes of the residual field, at rg = 7M, for puncture orders n — 2, 3 and 4. The upper plot 
shows that the small-m modes may take either sign. The lower plot shows power-law fall-off 5*^' ^ 0{m~'') 
at large m, with exponent C = 2 for the 2nd-order puncture, and C — 4 for the 3rd and 4th-order punctures. 
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FIG. 14: Modes of the radial (conservative) SF component, F™, at tq = 7M, for puncture orders n = 2, 3 
and 4. The upper plot shows that small-m modes may take either sign. The lower plot shows power-law 
fall-off F™ ^ 0{m~'^) at large m, with exponent C = 2 for the 2nd and 3rd-order punctures, and ^ = 4 for 
the 4th-order puncture. The dotted lines are reference lines oc and cx m~^. 
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FIG. 15: Modes of the angular (dissipative) SF component, i^™, at vq = 7M for puncture orders n = 2, 3 
and 4. The plot shows the modal contributions —F^ on a semi-log scale. It illustrates that (i) the modes 
f ™ do not depend on the order of the puncture to within numerical error (note that the points in the plot 
are superimposed upon one another), and (ii) the modal contributions diminish exponentially fast with m. 



converge exponentially fast in the latter case). For X"^ G {<I'^,F^}, we fit the simple power-law 
asymptotic model 



X™ = m-^ {A + B/m + C/m^ + ■■■), 



(105) 



where A,B,C, . . . are constant coefficients and C depends on the puncture order, as detailed at 



the start of Sec. IV B We are free to choose the number of terms N in the fit (105), and the 
part of the m-mode spectrum that we use for the fitting, i.e. mmin < m < rrimax, provided that 



N <m 



_ '"-max 



+ 1. Typical values in our analysis are N = 3, 



12, and m-n 



19. We 



split the sum into two parts. 



(106) 



m=0 



m=0 



m=r?lmax + l 



The first sum is found by adding the numerically-determined modal contributions. The second 



sum is found by analytically summing the fit formula, Eq. (105). The values of the fit parameters 
A, B,C, . . . depend somewhat on the set {N , mmin, nT'msLx} ■ By varying this set we may estimate 
the 'mode summation error'. Sample values of this error are quoted in the Tables presented in 
Sec. UVDl 



3. Mode cancellation error 



Figures 13 and 14 show that the modal contributions at small m may take either sign, and 
that the magnitude of individual modes can be substantially larger than the total sum. This is 
illustrated in Fig. 18, which compares the magnitude of the m = modal contribution (for 2nd and 



4th-order punctures) to the magnitude of the total mode sum. For both the field and radial SF, 
the m = mode is substantially larger in absolute value than the total. The total field (radial SF) 
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FIG. 16: Modal contributions to the radial component F^'^'^ for a range of orbital radii, with a 4th-order 
puncture. Plotted here on a log-log scale is the absolute magnitude of F™ as a function of m. The modal 
contributions diminish in magnitude as tq increases (see also Fig. 18), and it appears that all modes are 
scaled by approximately the same factor. For large to, the amplitude of the modal contributions falls off as 
~ TO~^. The asymptotic regime begins around m 10 for all radii. 



diminishes as Tq ^ (^o ^) ™ ^^e large-m limit (see ^24j), whereas the m = mode (and other modal 
contributions) are found to diminish far less rapidly. Hence the ratios $^^^/<I>/j and F^^^/F^'^^^ 
increase with tq. At ro = 30M, the former ratio is ~ 29 and the latter ratio is ~ 186. 

The phenomenon of cancellation between modes to leave a small remainder has a detrimental 
effect on the accuracy that can be achieved when computing the mode sum. At large rg we are 
reliant on delicate cancellations, with the total mode sum being orders of magnitude smaller than 
the m = contribution; this results in a relative error in the total mode sum that is much larger 
than the relative error in the individual modes. Unfortunately, this sets a practical limit upon 
the range of radii for which an accurate SF can be calculated using our version of the 4th-order 
puncture (Sec. IIIB4). It may be possible to find an alternate version of the 4th-order puncture 
which alleviates this problem somewhat; this issue is certainly worth further investigation. Note, 
however, that we do not anticipate mode cancellation to be a significant problem in the case of 
gravitational SF calculations, because in the gravitational case the radial SF diminishes as r^'^ at 
large ro (as opposed to in the scalar case). 



C. Computational resource 

The computational workload in the m-mode regularization scheme is 'embarrassingly parallel' 
in the sense that each run (i.e. each 2-l-lD evolution for given m, ro, {num.}) may be assigned to a 
separate thread, and little or no communication is required between threads. The mode sums are 
computed by post-processing the results from multiple runs. 

To run multiple threads in parallel, we made use of the IRIDIS 3 HPC resource. To obtain a 
SF estimate at a given radius, we typically compute 20 modes (m = 0, . . . , 19), at four different 
resolutions (e.g., n^es = 32, 48, 56, 64 with ares = 10). Thus we require approximately 80 nodes 
for each ro. For fixed grid dimensions, the runtime scales as ra^gg; hence the rires = 64 run takes 
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FIG. 17: Modal contributions to the angular component for a range of orbital radii. Plotted here on a 
semi-log scale is —F™ (as all modes are negative) as a function of m. The modal contributions F™ fall off 
exponentially with m and the decay rate increases with tq. 




FIG. 18: Comparing the m = mode ^^=° (left plot) and F™=° (right plot) with the totals and F^""^^, 
across a range of orbital radii tq. The dashed (red) line shows the magnitude of the monopole component 
for our 2nd-order puncture, and the dotted (blue) line shows the same for our 4th-order puncture. The solid 
(black) line shows the magnitude of the total radiative field (left) and radial SF (right), which scale as r^^ 
and (respectively) in the limit tq — oo. The ratios ^^^^/^r and F™^*'/^^"'^ increase as tq increases, 
and hence the accuracy of the mode sum degrades at large rg due to 'mode cancellation error' (see text for 
discussion). 



32 run. With U 



eight times longer than the n^ei 
approximately 12 hours. 

Additional resource is devoted to the m 



300M and 

^res — 64 a single run takes 



and m 



1 modes, to mitigate the problem of 
relaxation error (Sec. IV A 5). Typically we ran these modes up to t = lOOOM with a maximum 
resolution nres = 64 using the multigrid refinement scheme of Sec. IV A 5 c This illustrates a key 
flexibility of the m-mode scheme: the slow-decaying part of the initial junk is dominated by the 
lowest modes which can be handled separately. 
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Punc. order 


{M/qyF^°'' 




Tail contribution 


Large-m behaviour 


2nd 


7.86(7) X 10- 


-5 


31% 


0(m-2) 


3rd 


7.86(1) X 10- 


-5 


5.4% 


Cl(m-2) 


4th 


7.8507(3) X 10- 


-5 


0.2% 




f-domain 


7.850679 X 10" 


-5 







TABLE IV: Numerical results for the radial SF at rp = 7M . This table compares the results from imple- 
mentations of 2nd, 3rd and 4th-order puncture schemes against the frequency domain results (final row) of 
Diaz-Rivera et al. (see Table I in ^A\). The digit in paratheses indicates the estimated error in the final digit 
quoted; for example, 7.86(7) x 10"^ implies (7.86 ± 0.07) x 10"^. The third column ('Tail contribution') 
lists the proportion of the total radial SF which comes from the sum of the modes m > 15, i.e. the ratio 
Xlm=i6 Sm=o ^^^^ column indicates the asymptotic behaviour of modes at large m (see, 

e.g., Fig.fwl). 



D. Mode sum results 



In this section we present sample numerical results for the SF and radiative field, obtained using 
our (debut) implementation of the m-mode regularization scheme. Let us begin by considering 
results for the radial SF obtained with the 2nd, 3rd and 4th-order puncture schemes. Table IV 
shows numerical data for the radial SF F^^^^ for an orbital radius of rg = 7M, and compares with 
the frequency-domain /-mode results of |24| . As expected, the 4th-order scheme is most accurate 
and has the narrowest error bar, and the 2nd-order is least accurate and has the largest error bar. 
In the case of the 2nd-order puncture, the error is dominated by the 'tail-fitting error', i.e. the error 



in summing the large-m tail after fitting to an appropriate asymptotic model (see Sec. IV B 2 ). The 
tail- fitting error is large at 2nd-order for two reasons: (i) the tail decays slowly with m, as 0(m~^), 
and (ii) the contribution from the modes in the high-m tail represents a sizable proportion of the 



total (see Fig. 14). To illustrate the latter point, in Table IV we give the proportion of the total 
contained in the modes m > 15. We see that, although modes of the 3rd-order puncture also decay 
slowly, as 0(m~^), the tail at 3rd order has a smaller magnitude than the tail at 2nd order (~ 31% 
vs. ~ 5.4%), and hence the tail-fitting error in the 3rd-order result is commensurately smaller. In 
the case of the 4th-order puncture, the modes in the tail are rapidly-decaying [©(m"*^)], and the 



magnitude of the tail is small (see Fig. 14); hence tail- fitting error is much reduced and, in fact, 
no longer the dominant source of error. 

In Table |V] we present numerical results for the radiative field obtained via the 4th-order 
puncture scheme, for a range of orbital radii. The results shown in the second column were obtained 
by post-processing the results of multiple runs. We used unigrid runs up to tmax = 300M for modes 
m = 2, . . . , 19 at a range of resolutions n^es = 32, 48, 56, and 64 with Ores = 10. We used the 



data (t 



IV A 5 


), we 1 


Sec. 


IVA5c) 



various resolutions to extrapolate to zero grid spacing (see Sec. IV A 2). To mitigate the relaxation 

lOOOM on a three-level 
fitted the late-time 



and m 



1 up to imaj 
= 64, and for m 



900M to lOOOM) with the appropriate power-law relaxation model. To sum the large-m 

19 with a three-term 



12,... 



tail (Sec. IV B 2), i.e. the modes m > 19, we fitted the modes m 
model Am~^ + Bm^^ + (jm-6_ Estimates of the residual errors that remain after performing these 
steps are given in the final three columns. We find that, although the residual errors are broadly 
similar in magnitude, the residual relaxation error remains the largest. 

In Tables |VI| and |VII we present numerical results for the conservative and dissipative com- 
ponents (respectively) of the SF, for a range of radii. The radial component of the SF, shown 
in Table |VI[ was computed in a similar manner to ^r, that is, by post-processing the results of 



multiple runs to minimize discretization error ( IV A 2 ) , relaxation error ( IV A 5 ) and tail- fitting 
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error (IV B 2). The angular component of the SF, shown in Table VII, was simpler to compute 
because in this case it was not necessary to model the large-m tail, since the modal contributions 



decay exponentially- fast (see Fig. 15). In addition, relaxation error is less significant for F^^^ 
because the slowly-relaxing m = mode of F^^^ is zero. We remind that the temporal component 
of the SF, Ff^^\ may be found directly from F^^^ using (38). 

Tables V VII demonstrate that the m-mode methodcan yield highly-accurate SF estimates 
for circular orbits in the strong field (i.e. 6M < ro < lOM). For example, at tq = 6M, the 
conservative part of SF (i.e. F^'^^^) is in error by one part in ~ 5 x 10^, and the dissipative part 
of SF (i.e. F^"^^, Fi^^^) by one part in ~ 5 x 10^. This level of accuracy would not have been 
possible without careful modelling and mitigation of the sources of error described in the previous 
sections. Unsurprisingly, the results for F^^^, whose large-m modes converge exponentially fast, 
are substantially more accurate than the results for F^"^^^, whose large-m modes exhibit power-law 
decay, 0{m~'^). 

It is clear from Tables V -VII that the results of the m-mode method degrade in accuracy as tq 
increases. For example, the relative error in our result for F^'^^^ at tq = 30M is approximately 4500 
times greater than the relative error at tq = 6M. A loss of relative accuracy with radius can be 
explained from two points of view. Firstly, the magnitude of the singular part of the field (and thus 
the magnitude of the retarded field outside the worldtube, which is related to the magnitude of 
the discretization error) depends only weakly upon the orbital radius. In other words, the absolute 
error in our simulations depends only weakly upon tq, and therefore, since the radial SF falls off 
very rapidly (as r^^)^ the relative error in our results increases rapidly with tq. An alternative 
point of view is to see the relative loss of accuracy as a consequence of 'mode cancellation error', 
which was described in Sec. IVB 3[ Fig. |18| shows that the ratio of the magnitude of a typical mode 
to the total mode sum increases rapidly with rg; hence small relative errors in individual modes 
may become large relative errors in the total. 

The loss of accuracy at large tq demonstrated here is not a particular concern to the prospect 
of accurate time-domain SF calculations, for two reasons. Firstly, in the more interesting case of 



the gravitational SF, the radial component of the SF falls off only as 



rather than 



and 



so we expect the loss of accuracy to be far less pronounced in that case. Secondly, we note that 
complementary approaches, such as Post-Newtonian methods, are well-suited to modelling orbits 
in the weak-field (large-ro) regime, whereas the primary goal of the SF approach is to accurately 
describe EMRI physics in the strong-field (small-ro) regime. 



V. DISCUSSION AND CONCLUSIONS 



In the foregoing sections, we have presented details of the first implementation of the m-mode 
regularization scheme for SF calculations. Our objective has been to establish the scheme, first 
proposed in [591 160], as a practical alternative to (i) /-mode regularization in cases where separation 
of variables seems infeasible, e.g., gravitational SF on Kerr in the Lorenz gauge, and (ii) 3 -|- ID 
evolution schemes (see, e.g., [31]). We have demonstrated here that (at least in the simple case 
of circular orbits on Schwarzschild) the results from /-mode and m-mode regularization are fully 
consistent, to within the limits of numerical accuracy. This work provides reassuring evidence that 
the m-mode scheme is well-founded and has been correctly implemented. 



To attain accurate results (such as those presented in Sec. IV), we were compelled to improve 



our understanding of a number of issues, including (i) the effect of puncture order upon the rate 



of convergence of the m-mode sum (Sec. II G and IVB), (ii) the stability and convergence of the 



finite difference scheme (Sec. IIIC5), (iii) the influence of various sources of error upon numerical 



accuracy (Sec. HIE and IV A), and (iv) the use of a simple mesh refinement algorithm to improve 
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TABLE V: Numerical results for the field <i>fl for a range of orbital radii. The second column gives the 
numerical result from our implementation of the time-domain m-mode method using the 4th-order puncture, 
a maximum resolution n^cs = 64, arcs = 10, unigrid runs with tmax = 300Af for modes m = 0, . . . , 19, and 
multigrid runs for modes m — and m — 1 with tmax — lOOOA/. The third column quotes the (highly- 
accurate) frequency-domain /-mode method results of [24 for comparison. The remaining columns give 
estimates of sources of numerical error. The relaxation error ('Relaxation') was estimated from alternative 
extrapolations to infinite time for the to = mode, using riics = 64 multigrid data up to t = lOOOM. The 
discretization error ('Discret.') was estimated by summing in quadrature the modal discretization errors 
found by comparing alternative extrapolations to infinite resolution based on data at n^cs = 64,56,48,32. 
The large to tail-fitting error ('Tail fit') is an estimate of error in summing modes in the large-TO tail above 
TO > 19 by using a fitting model Am^^ 4- Bm~^ + Cm~^ . These estimates suggest that relaxation is the 
dominant source of error in for all radii. The error bar on the final result (parenthetical figures in the 
second column) was found by combining these error estimates in quadrature. 
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5.2x10" 


-4 


14 


2.721(3) 


X lO-'' 


2.72008 X 10" 


-6 


2.5x10-4 


2.3x10" 


-4 


9.6x10" 


-4 


20 


4.96(3) 


X 10-^ 


4.93790 X 10" 


-7 


5.5x10-3 


2.7x10" 


-3 


2.9x10" 


-3 


30 


8.5(8) 


X 10-« 


7.17192 X 10" 


-8 


2.0x10-2 


7.9x10" 


-2 


6.5x10" 


-2 



TABLE VI: Numerical results for the radial SF for a range of orbital radii. The second column gives 
the numerical result from our implementation of the time-domain TO-mode method. The structure of the 
table is similar to that of Table |V] and the same numerical parameters have been used. 













ro/M 


TO-mode, t-domain 


/-mode, f-domain [24] 


6 


-5.30423170(3) 


X 


10-3 


-5.30423170 x 10" 


-3 


7 


-3.2731229(1) 


X 


10-3 


-3.27312280 x 10" 


-3 


8 


-2.2111614(1) 


X 


10-3 


-2.21116134 X 10" 


-3 


10 


-1.1859260(2) 


X 


10-3 


-1.18592599 x 10" 


-3 


14 


-4.838491(2) 


X 


10-4 


-4.83849328 x 10" 


-4 


20 


-1.92445(2) 


X 


10-4 


-1.92444253 x 10" 


-4 


30 


-6.8629(3) 


X 


10-5 


-6.86315934 x 10" 


-5 



TABLE VII: Numerical results for the angular component of SF, \ for a range of orbital radii. This 
data were obtained with the same numerical parameters as in Tables |V] and |VI[ The third column shows 
a comparison with the results of |24) . The error here is predominantly discretization error, arising from 
extrapolation to zero grid resolution. 
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computational efficiency (Sec. IV A 5c). A crucial step forward was taken in moving from a 2nd- 
order puncture scheme (described in [60j) to a 4tli-order puncture scheme. We believe that the 
4th-order scheme will be the foundation for a range of m-mode implementations. 

An obvious way to achieve greater accuracy is to increase the grid resolution. However, given 
that the scaling of runtime with resolution is problematic (runtime oc nj?gs), more advanced methods 
may yet be needed. One possibility is to use a higher-order (e.g., 4th-order) finite-differencing 
scheme. A subtlety here is that global 4th-order convergence may be difficult to achieve, given 
the limited differentiability of the effective source on the worldline. Case-by-case treatment of 
finite difference molecules near the worldline will be necessary, but somewhat arduous. Another 
possibility is to employ a systematic adaptive mesh refinement scheme a la Thornburg [62]. This 
is a promising route for the future, but one which will require a much more sophisticated code 
architecture. 

It is possible that the use of higher-order punctures may become feasible in future. Higher- 
order punctures would improve the power-law convergence of the m-mode tail of the conservative 
component of the SF, and generate an effective source which is smoother and flatter near the 
worldline. However, a naive implementation would face at least two difficulties: (i) computing 
the d'Alembertian of the puncture algebraically leads to an extremely cumbersome expression at 
high orders, and (ii) computing a finite source from a divergent puncture may rely on delicate 
cancellations between terms, a problem which (without care) will get worse at higher order (but 
see Ref. [61] for possible resolutions of these problems). For practical reasons, we believe the 
4th-order puncture may represent a 'sweet spot' since it is the lowest-order puncture to exhibit 
m~'* convergence for the SF modes, which renders the tail-fitting error sub-dominant. To obtain 
improved convergence for the SF modes, one must go up to sixth order. 

An intriguing possibility is that calculations of 'infinite order' (i.e., quasi-exact) punctures may 
become possible by, e.g., obtaining U and V in the Hadamard form (11) by integrating transport 
equations along the family of geodesies that join a field point to the worldline [7S]. In this case the 
effective source would be identically zero, and all components of the SF would exhibit exponential 
convergence with m. See Ref. [61] for further discussion of this idea. 

In this work, we have solved the wave equation in the region exterior to the horizon, using a 
tortoise coordinate and a uv grid. This approach has the benefit of simplicity, in that boundary 
conditions are not required at the horizon or at spatial infinity. However, the uv method has several 
drawbacks: (i) doubling the simulation time quadruples the grid area and thus the run-time, (ii) 
much runtime is 'wasted' in the near- horizon regime of small- r*, (iii) the reflection of 'junk radiation' 
from large r leads to rather slow power law relaxation. An alternative spacetime slicing, such as 
the asymptotically-null hyperboloidal slicing proposed in [76j, may bring advantages. The benefits 
of such a slicing for SF calculations were recently highlighted [77]. This example of 'technology 
transfer' highlights the fact that many of the techniques routinely employed by numerical relativists 
could also benefit time-domain SF simulations. 

A companion paper (in progress) will describe the first implementation of the m-mode scheme 
for scalar-field SF on circular orbits of the Kerr spacetime. As the m-mode scheme is axisymmetric 
by construction, the method outlined here actually requires little modification. In the forthcoming 
work, we will focus only on the additional issues that arise, such as the problem of finding a 
stable finite-differencing method. Fortunately, we may test the accuracy of our implementation by 
comparing against recent /-mode results j56] . 

In progressing the m-mode scheme, we remain mindful of three open challenges for the future. 
The first challenge is to compute the gravitational SF in the Lorenz gauge. We hope to break 
new ground by examining the gravitational SF for circular orbits on Kerr within the m-mode 
scheme. Such a calculation has not been possible with the Z-mode scheme, due to the apparent 
inseparability of the field equations in the Lorenz gauge. Of course, it may be possible to work in an 
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alternative gauge; recent progress on a radiation gauge calculation is described in [35]. The second 
challenge is to adapt the scheme to treat highly eccentric or unbound orbits, which are beyond the 
scope of frequency-domain approaches. The third challenge is to achieve accurate, self-consistent 
long-term orbital evolutions of EMRIs using gravitational SF calculations. We refer the reader to 
Refs. [TTl [50l [78] for first steps in this direction. 
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Appendix A: Toy model for the residual field: Details 



Here we show that the function H^^^ introduced in Sec. II G Eq. (39), has an m-mode contri- 
bution given by Eq. (41). 

Although it is non-smooth at (/? = 0, -ff'"' is differentiable an infinite number of times in a 
piecewise sense. Hence we may express its derivatives in terms of distributions. In particular, near 
= 0, its (n + 2)th derivative has the expansion 

^N(n+2) ^ 2hn6"i^) + 2hn+l6'{ip) + 2K+26{ip) + sign(99)(. • • ). (Al) 

Here 6{-) is the Dirac delta distribution, a prime denotes differentiation with respect to ip, and (• • • ) 
represents a regular Taylor expansion in ip about ip = 0. We have made use of the distributional 
identities ip/\(p\ = 2Q{<p) - 1, Q'{(p) = 6{(p), ^^5^^\ip) = k\{-lf6{<p) and >p^5{ip) = for integer 
A; > 0. Note that the globally defined function 

^^[n](n+2) ^ jj[n](n+2) _ 3 + hn+l5' {^) + /i„+2<5(^)] (A2) 

is hounded in magnitude everywhere on — vr < < vr. 

Now consider the m mode coefficient of H^'^\ given (for m 7^ 0) by 

2ttJ_, ^ 27r(im)"+2 ^ 

1 ^^[„](n+2)g-*m^^^ 



1 

7r{im) 



+ ^,,^,n+2 / [hnS"{v) + K+l5\^) + K+25{^)\ C'^^^d^, (A3) 



where in the second equality we integrated by parts n + 2 times (note that we used the condition of 
smoothness across ip = — vr, vr to eliminate boundary terms). Consider the first integral in the final 
expression: Since ["l("+2)g-«m'/'| bounded on — vr < (p < tt, the magnitude of this integral can 
be bounded by Cm~^~'^ with some positive (m-independent) constant C. The second integral in 
the final expression is readily evaluated in explicit form. Altogether we get (for m 7^ 0) 

(zm)" (2m)"+i (zm)"+2 ^ ' 



(zm)" (im)"+^ 
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where in the first hne ©(m"""^) represents the contribution from the (5ii'N("+2) integral in 
Eq. (A3), and in the second line we have absorbed the oc hn+2 term within ©(m"""^). Inserting 
(A4) into the modal contribution formula (33) leads directly to Eq. (41). 



Appendix B: to- mode decomposition in the 2nd-order scheme 

The individual modes of the puncture field are obtained via integrals, 

1 



$^ = — / $pe-*™^d(^. (Bl) 



— 7r 



With <I>-p as defined in (52), and after the replacement in Eq. (55), these integrals can be performed 



analytically, in a similar manner to [SS]. We find 

= Wt^^ [p^^(p)ellipK(7) +p^(p)ellipE(7)] (B2) 



where 

and the quantities A and B are simply 

A = Prr5r'^ + Pee6e'^ + Qrr5r^ + Qee5r5e'^, B = P^^ + Q^^6r, (B4) 



with coefficients Pij and Qij defined in Eq. (51) and (54). Here ellipK(-) and ellipE(-) are complete 



elliptic integrals of the first and second kinds, respectively, defined by 

ellipK(ifc) = / sin^ x)'^^^ dx, ellipE(fc) = / (l - A;^ sin^ x) ^''^ dx. (B5) 

Jo Jo 

The polynomials p^{p) and p^{p) were given explicitly for m = 0, . . . , 5 in the tables of Appendix 
A of [59J (note that our p plays the role of p in Ref . |59] ) . Polynomials for m > 5 are straightforward 
to calculate using a symbolic algebra package. 

The m modes of the effective source are obtained via the integral 

STs =^ r Sefie-^'^'^dyp, (B6) 



where S'efj = S — □'^-p. These integrals may be expressed analytically in the form 



where 



ztt 



51 = (r - M)r-^X{r, B) + f{r) {P„ + ZQrM + r'^ (1 + 89 cot 9) {Pee + Qee^r) , (B8) 

52 = {r-^sm-^9-ujyfir))B-2ir-M)r-^Q^^, (B9) 
^3 = -3f{r){Xir,9)/2-Q^^f -3r-^69^Pee + Qe9Sr)\ (BIO) 
54 = -W^ (r-2 sin-2 _ oo^fir)) + 3/(r)Q2^, (Bll) 

^5 = -3fir){X{r,9)/2 + Q^^f -3r-^69^Pee + QeeSr)\ (B12) 
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and 



(B13) 



In Eq. (B7) the quantities are 

ep3 = ^ [pri^(p)ellipK(7) + p-^pTj,{p)empEij)] , (B14) 

cos6ipe-'^'^di6v) = [p^;,(p)ellipK(7) + p'^p- (p)ellipE(7)] , (B15) 



jm 
-'1 



jm 
-'2 



6ip 



7 



^2^5/2 



[pS^(p)ellipK(7) + (p)ellipE(7)|pi6) 



{dip) e-'^^^^didip) 



7 



^5/2 



[prx(p)enipK(7) + p-2pS;(/3)ellipE(7)] , (B17) 



9^) 



7 



^5/2 



fe(p)ellipK(7) + p-2pri?(p)enipE(7)] . (B18) 



Note that I™, and I™ are the same integrals as defined in Eqs. (46), (49) of [59] but with P^^ 
replaced by B and with replaced by A/{4:B). Note also that we have introduced new definitions 
for Is and I5. The polynomials P2K' P2Ey pTk p^E were given for m = 0, . . . , 5 in 

Appendix A of [59]. The polynomials and again for m = 0, . . . , 5, are given in 



Table VIII here. Polynomials for m > 5 were calculated using a symbolic algebra package. 



Appendix C: Stability of the finite difference method 



As discussed in Sec. Ill C 5 in vacuum simulations we observed a numerical instability arising 
first near the poles, with a short wavelength 2A in the 9 direction and an exponentially-growing 
amplitude. The origin of the instability can be better understood by applying a von Neumann 
stability analysis (see, e.g., [71]) to the finite difference equations (87) in vacuum (i.e., for = 0). 
Let us consider some numerical 'noise' on a timeslice t = tc with a short angular wavelength of 
2tt/h (with K S> 1) and an amplitude e^,- If it turns out that the finite difference method amplifies 
this noise exponentially, then we expect the method to be unstable. Let us begin with an ansatz 



*"^(ti,r,-,0fc) = e,e 



(U-t,)/(h/2) 



exp {inOk) 



(CI) 



where ti , rj , and Oj, are the values of coordinates at points 1-8 in the grid cell shown in Fig. [3| Note 
that here we have ignored variation in the r direction to focus only upon the angular instability. 



Next we insert (CI) into the finite difference equation (87) to obtain a quadratic equation for the 



(complex) amplification factor ^, 



+ re + 1 = 0, 



where 



T 



-2 + 



4^2 ^2 



2(1 - cos(kA)) - iA cot dk sin(KA) + A^ 



2M 



+ 



m 



sin^ Ok 



(C2) 



(C3) 



The general case is difficult to analyse, so let us focus on the relevant case of the short- wavelength 
angular mode with k = vr/A, in the grid cell closest to the North pole with 0^ = A. This is precisely 



43 



m pf^ip) 
-21 

2 ^(6V + 40p2-l) 

3 ^ (512p6 + 640/ + 180p2 - i) 

4 (I6384p** + 29696p6 + 16512^-1 + 2720/92 - 5) 

5 (655360p^° + 1540096p^ + 1263104p6 + 418688p^ + 45500^^ - 35) 



TO p^Eip) 

1 -i(4p4 + p2_2) 

2 -^(64pfi + 72^-1 + 7^2 -2) 

3 (512p« + 896p6 + 404p4 + 17 - 2) 

4 - fl6384pio + 37888p8 + 28288/9^ + 6944^^ + 155^2 _ ^q) 

5 -8io (655360p^2 ^ i867776p^° + 1910272p^ + 822912^^ + I26876p^ + 1715^^ - 70) 



m p^^f (p) 

1 -^(4p2+5) 

2 (64p4 + 88p2 + 23) 

3 - ^ (512p6 + 896p4 + 436p2 ^ 53-) 

4 - j|g h6384p« + 35840p6 + 25728^'' + 6752p2 + 475) 

5 -sin (655360p^° + 1736704^^ + 1656320p'^ + 683648/?^ + 113852^2 + 5215) 



24 
J_ 
24 
1 



4p4 



-1) 

7p' 



1 



. 24 v6V + 120p* + 55p2 + 1) 

3 ^ (512p8 + 1152p6 + 788p4 + 151p2 + 1) 

4 i fl6384pio + 44032/9* + 40576p'' + 14432^* 
^ ^ (655360pi2 + 2064384pio + 2401792^* + 1247616^^^ 



1499p2 



840 



h5) 

272412p4 



17669/92 + 35) 



TABLE VIII: The polynomials p™^ and p™^ appearing in Eqs. (BI6I and (BI8I, for m = 0, . . . , 5 



the mode which was observed to be the source of the instability in our numerical implementation. 
In this case, T is the real quantity 



T 



-2 + 



4^2 ^2 



(C4) 



If —2 < T < 2 then the roots ^ of (C2) are a complex-conjugate pair with |^| = 1 and so we expect 



the method to be stable. Conversely, if |T| > 2 then at least one root of (C2) has a magnitude 
larger then unity, and we expect the method to be unstable. The requirement |T| < 2 in (C4) 



leads immediately to the stability condition (89). 
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